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THEORY OF DISORDERED HEISENBERG FERROMAGNETS* 
by Robert M. Stubbs 
Lewis Research Center 

SUMMARY 

The theory of Heisenberg ferromagnets with distributed disorder is developed by a 
method using double-time, temperature-dependent Green's functions. Disorder is in- 
troduced by allowing the exchange interactions between spins to deviate randomly from 
the mean interaction. The disorder is characterized by a parameter p which is pro- 
portional to the mean- square deviation of the exchange interactions from the mean inter- 
action. The equations of motion for Callen-type Green's functions are solved by using 
Tyablikov decoupling, and an ensemble average is performed over systems with similar 
disorder to provide an ensemble-averaged Green's function. From these disordered 
Green's functions the densities of spin wave states are derived which, in turn, are used 
to calculate the magnetic properties of disordered systems. Several specific systems 
are investigated including disordered simple cubic, body- centered cubic, and face- 
centered cubic systems with various spin values. The theory is also applied to mixtures 
of the chalcogenides of europium, which are good examples of Heisenberg ferromagnets. 
Disorder is shown to have a marked effect on the density of spin wave states. Modest 
values of the disorder parameter can produce relatively large changes in the state den- 
sity - in the form of enhancement of the low-energy densities and extension of the energy 
band to higher values. The increase of the density of low-energy states due to disorder 
is of the order (1 - p)~^^ , leading to a corresponding increase in the low- temperature 
specific heat. The spontaneous magnetization of a disordered ferromagnet decreases 
with rising temperature more quickly than for a crystal, and the Curie temperature is 
shown to decrease linearly with disorder as 1 - p. Calculations for the low- 
temperature region and for the higher temperature paramagnetic phase show that dis- 
order effects are more pronounced in the low- temperature, ferromagnetic phase. 


*The material in this report was submitted as a thesis in partial fulfillment of the 
requirements for the degree Doctor of Philosophy at the University of Toledo, Toledo, 
Ohio, May 1972. 



INTRODUCTION 


Until recent years, solid-state physics has been a narrower field of study than its 
name suggests. Perhaps crystal physics would have described the field more accurately 
since nearly all published literature on the physics of solids, both theoretical and ex- 
perimental,. was confined to systems with translational symmetry. This is not surpris- 
ing since the mathematical simplifications inherent in lattice symmetries are consider- 
able. Nor did this state of affairs present much of a limitation since most materials of 
technological importance had a crystalline structure. 

In 1965, the first amorphous ferromagnet was discovered (ref. 1) and several more 
were identified soon afterward (refs. 2 to 5). When they arrived on the scene, there 
was virtually no theoretical understanding of what to expect when disorder was present 
in a magnetic system. The fact that ferromagnetism is a cooperative phenomenon 
makes even more intriguing the question whether disorder would significantly affect a 
ferromagnet and, if so, in what manner. In a sense the ions in a ferromagnet are more 
strongly influenced by the other particles in the system than in an uncooperative phase 
such as a paramagnet. In a ferromagnet the magnetic moments act in concert to pro- 
duce a more or less parallel alinement, while the moment directions in a paramagnet 
are random and the net magnetic moment vanishes. One might ask if the introduction 
of disorder would have drastic results on such cooperative effects. In its present state, 
the theory of magnetism is not sufficiently developed to give a definitive answer to these 
questions. It is the intent of this report to rectify some of this incompleteness by ex- 
panding and generalizing the theory of mangetism to include the effects of disorder. 

In this report the problem of disorder in a ferromagnet is treated by a Green's 
function theory that incorporates disorder by allowing a randomness in the exchange 
interactions coupling the magnetic moments. The advantage of the Green's function 
method is that it gives good results over the entire temperature range (ref. 6). In the 
interesting regions near T = 0 and near T = T^,, the Curie temperature, and in the 
high- temperature paramagnetic region. Green's functions for disordered systems are 
derived and used to calculate such properties as the density of spin wave states, spon- 
taneous magnetization, Curie temperature, and specific heat and their dependence on 
disorder. These calculations are applied to systems of various values of spin and vari- 
ous structures. 

Some background material for this study is provided in the form of a brief sketch of 
the present state of the art. It is of interest to know which types of disorder have re- 
ceived attention, what techniques have been used, and which problems have been solved 
and which have not. The formal theory of double-time, temperature- dependent Green's 
functions is presented. These are the Green's functions most useful in statistical phys- 
ics. These Green's functions are then used to treat the problem of a disordered Heisen- 
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berg ferromagnet. Specific disordered systems are studied; and, finally, general con- 
clusions are drawn from the various systems investigated. The appendixes contain 
some of the more lengthy derivations and calculations. 


EARLIER WORK 

The modern theory of magnetism has developed rapidly since 1926 when Dirac 
(ref. 7) and Heisenberg (ref. 8) independently discovered the concept of exchange. 
Classical physics had no way to explain the strong interactions between spins needed to 
account for ferromagnetism. It was shown that exchange was a purely quantum mechan- 
ical effect with no classical analog and depended on the overlap of the wave functions of 
the electrons whose spin accounted for the magnetic moments. The energy of interac- 
tion between a spin at site i and a spin at site j was shown to be . - J(ij)Si • Sj, the so- 
called Heisenberg energy, where J(ij) is the exchange integral. Depending on how com- 
plicated a ferromagnetic system is considered, the Hamiltonian might contain several 
more terms. But the Heisenberg term is basic and can be derived from the fundamental 
Coulomb interactions of the electrons together with the Pauli exclusion principle. 

Bloch (refs. 9 and 10) analyzed the Heisenberg model of a ferromagnet by using the 
concept of spin waves. A spin wave consists of a single reversed spin distributed co- 
herently over a large number of otherwise alined atomic spins in a crystal. Dyson 
(refs. 11 and 12) considered the interactions between these spin waves and derived a 
low- temperature series expansion for the magnetization in powers of the absolute tem- 
perature T. Other temperature regions were investigated by alternate methods. For 
temperatures near the Curie temperature T^, molecular field theory (ref. 13) was 
used; and for high temperatures, perturbation theory (refs. 14 and 15) was employed to 
give an expansion in 1/T. In 1959, Tyablikov (refs. 16 and 17) showed that double- time, 
temperature- dependent Green's functions could be used to describe the Heisenberg fer- 
romagnet over the entire temperature range. His results for a spin -1/2 system showed 
agreement in the main terms with all three methods. The Green's function theory was 
extended to treat systems of spin greater than 1/2 by Tahir-Kheli and ter Haar (ref. 18) 
and by Callen (ref. 19). Several attempts (refs. 20 and 21) to refine their work have 
been made since, most having to do with improved decoupling procedures. 

The first investigation of a noncrystalline ferromagnet to appear in the literature, 
either theoretical or experimental, was a 1960 paper by Gubanov (ref. 22) in which he 
used a semiclassical method to demonstrate the possibility of th.e existence of amor- 
phous ferromagnets . Since that time, there has been a steady growth in the study of 
various types of disorder in magnetic systems. One of the first problems treated was 
the case of an ordered Heisenberg ferromagnet in which one spin is replaced by an 
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impurity spin differing from the host atoms in either spin magnitude or exchange coup- 
ling. This system has been examined in detail by Wolfrom and Callaway (ref. 23) and 
by several others (refs. 24 and 25), all of whom showed the existence of a spin wave 
mode localized on the impurity, in addition to a modification of the spin wave band of 
the host lattice. Hone and Vogelsang (ref. 26) demonstrated that weakly coupled impuri- 
ties produce low-lying spin wave resonances leading to specific- heat anomalies at low 
temperatures. These results are valid only for low concentrations of substitutional im- 
purities since interactions between impurities are neglected. 

Handrich (ref. 27) studied amorphous Heisenberg ferromagnets by using molecular 
field theory, which predicted a decreased spontaneous magnetization arising from fluc- 
tuations in the structure of the system . A shortcoming of this approximation is its fail- 
ure to show any dependence of the Curie temperature on these structure fluctuations. 
Later work (refs. 28 and 29) within the same molecular field approximation using com- 
puter experiments showed increases in the Curie temperature brought on by randomness 
in the exchange interactions. This behavior is not that to be expected from other stud- 
ies (refs. 30 and 31) that show disappearance of ferromagnetism above critical concen- 
trations of defects. 

Recently, Montgomery, Krugler, and Stubbs (ref. 32) treated distributed disorder 
in a spin -1/2 Heisenberg ferromagnet by a Green's function method. This was the first 
quantitative demonstration that disorder decreased the Curie temperature. Because 
that work represents the preliminary part of this study, further comments on it will be 
saved for later sections of this report. 


GENERAL THEORY OF TEMPERATURE-DEPENDENT GREEN’S FUNCTIONS 

Often one wishes to calculate the expectation value of products of operators, and a 
powerful and elegant technique that accomplishes this involves the use of Green's func- 
tions. We will consider the two kinds of Green's functions which are most convenient 
in the statistical treatment of magnetic systems: the retarded Green's function G r (t, t '), 
and the advanced Green's function G„(t,t’). The definitions (ref. 6) of these and the 
notation used in their regard are 

G r (t,t’)s«A(t);B(t')» r =-i0(t- t')<[A(t),B(t’)]> (1) 

G a (t,t') = <<A(t);B(t'))) a = i0(t’ - t)<[A(t),B(t')]> (2) 

where 6 is a step function, 
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In general, the operators A(t) and B(t) are particle creation and annihilation operators 
or products of these. 

There are some properties of these Green's functions that are worthy to note before 
proceeding. First, neither of them is defined when t = t’ because of the discontinuous 
function 0. Secondly, G„(t, t’) and G_(t,t') are dependent on time only through t - t ' . 

r d 

This is seen by writing the definitions in explicit form, using equations (4) and (5) and 
making use of the commutability of operators under the trace sign. And, finally, the 
averages involved in equations (1) and (2) are not taken over the vacuum state of the 
system but over the grand canonical ensemble. And it is here that the Green's functions 
receive their temperature dependence. 

The advantage of Green's functions is that, once their solution is found, they can be 
used to calculate various correlation functions, from which all the properties of interest 
in a system can be derived. The equations relating correlation functions and Green's 
functions are developed in a later discussion of spectral representations. Now, we are 
interested in writing the equation of motion for the Green's function and finding its solu- 
tion. 

The time rate of change of the retarded Green's function is 
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(V 


iK — G (t, t') = K — (t- t')<[A(t),B(t’)]> - i0(t- t')<[[A(t),3C],B(t’)]> 
dt r dt 


Noting that 


A 0 (-t) = - A e( t ) (8) 

dt dt 

we see that both the retarded and advanced Green’s functions obey the same equation of 
motion and the subscripts can be dropped. An integral representation for 6 involving 
the Dirac delta function is 



( 9 ) 


allowing equation (7) to be rewritten as 

ifi— ( <A(t); B(t’)» = K6(t- t’)< [A(t),B(t’)]> + <( [A(t),3C]; B(t’)>> (10) 

dt 

This equation contains two Green’s functions plus an inhomogeneous term involving 
a delta function which recreates the form of the usual classical Green's function equa- 
tions. Classically, Green's functions are used to obtain the field caused by a point 
source or a distribution of point sources, and the solution is generally an integral repre- 
sentation involving the Green's function. In these cases the Green's function is a solu- 
tion to a differential equation having an inhomogeneous delta function term. It is the 
similarity of these equations with those like equation (10) that brought about the label 
"Green's function" for the functions used in quantum field theory. 

However, the Green's function on the right side of equation (10) has, in general, 
more terms than the original Green's function; that is, the operator [A(t),3€] has a high- 
er number of terms than does the operator A(t). To solve equation (10) exactly, one 
needs to know this higher Green's function, (( [A(t),3€]; B(t')}) , which requires the solu- 
tion of the equation 

ifi~ (< [A(t),3C]; B(t'))> = B5(t - t')< [[A(t),3C], B(t’)]) + « [[A(t),3C],3C]; B(t')>> (11) 

dt 

This introduces still a higher Green’s function, (( [[A(t),3C],3C]; B(t ))), and it is clear 
that we can generate an infinite set of coupled equations involving a hierarchy of Green's 
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functions of the form (< [. . . [[[A(t),3C],3C],3C]. . . ];B(t’))). Although the exact solution 
for ((A(t); B(t'))) requires solving the chain of equations, one can sometimes make an 
approximation that decouples the chain and leaves a finite number of equations which can 
be solved. This is, in fact, the procedure we employ later when we discuss decoupling 
more fully For the present, we assume that the chain typified by equations (10) and 
(11) can be decoupled and a solution found for the Green’s function. The question now is 
how can G r (t,t') and G a (t,t') be used to find the quantities of interest in the system un- 
der consideration. To calculate these properties from the Green's functions, use is 
made of correlation functions and spectral representations which we discuss in some 
detail in the following paragraphs. 

When written out explicitly, the equations for Gft,t') and G„ (t, t ’) contain terms 
like (A(t)B(t')) , which are averages over the grand canonical ensemble of products of 
Heisenberg operators. They are called correlation functions and are extremely impor- 
tant in statistical physics. When t 4 t', these averages are the time correlation func- 
tions which are useful in the calculation of transport properties. At equilibrium, the 
time correlation functions depend only on t - t' as do the double-time Green's functions, 
again because of the commutability of operators under the trace sign. When t = t', 
then, 


<A(t)B(t)) = (A(O)B(O)) = (AB) (12) 

These are the more usual correlation functions which are used to evaluate the average 
value of the dynamical quantities associated with the operators. 

To develop the spectral representations for the correlation functions, we carry out 
explicitly the operation implied by the angular brackets, ( . . .) , using the sum over 
eigenstates of the Hamiltonian of the system. We assume the set of states | / u) , where 

3C|m> =E j Jm> (13) 

to be complete and orthonormal so that we can rewrite equation (4) as 

-E /kT 

(A(t)B(t’)) =Q" 1 £(M|A(t)B(t')|ju)e M (14) 

P 


where Q is the partition function, 


Q , Tr(e- 5C / kT ) 


( 15 ) 


The Dirac notation on the right side of equation (14) should cause no confusion with the 
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angular averaging brackets on the left. We shall expand equation (14) further in the 

same basis, T"] f v) ( v [ , and use relation (5). 
v 


-II 


(A(t)B(t')) = Q" 1 y < M | A(0) 1 < i/ 1 B(0) | M> e 

M v 


-i(E M -E p )(t*-t)/R -E^/kT 


(16) 


Similarly, 


( B(t’)A(t)> = Q 


■II 


-i(E -E )(t-t')/h -E AT 
</x|B(0)|v)<i/|A(0)|M)e M e M 


l± v 


(17) 


By defining a spectral density, 


J AbM * Q- 


-II 


( p | A(0) | v) ( v | B(0) | \x) e 


E > T a <„, 


h h 


p V 


we can rewrite equations (16) and (17) more compactly as 


✓*■30 

(A(t)B(t')) = / I AB (w)eRW,/kT dw 

«/_ oo 


(19) 


<B(t’)A(t)> = / I AR ( W )e _iw(t_t,) dw 

- 00 


AB V 


( 20 ) 


Expressed in this form, the spectral density I^g(^) is seen to be the Fourier transform 
of the time correlation function. 

The spectral representations of the time correlation functions, equations (19) and 
(20), can be used to write the spectral representations of the retarded and advanced 
Green's functions. First, we define the time Fourier transform of the Green's functions 
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cyt 


- t') = f ftW(=)e' i(E/il)(t " t,) dE 
L W a 


( 21 ) 


if 


G (E) = -?- / G_(t - t')e^ E//R)t dt 
r 2~ 1 r 


More explicitly, the time transform of the retarded Green's function is 


G r (E) = — f -i0(t)|( A(t)B(t')) - ( B(t')A(t))J e 1 

271 Jo O 


(E/R)t 


dt 


or, by equations (19) and (20) 


( 22 ) 


,(E) = — f dt 9(t) f dcjI AB (o))e i S E/K)_w ] t (e na;//kT - 1) (23) 

2lt Joo J-o o 


At this point we notice that the step function 0(t) can be represented in integral form 


, f -icu’t 

0(t) = lim — / dec' 

e— 0 27 t I u)' + ie 

c/- OO 


e > 0 


(24) 


That this integral representation has the properties of the discontinuous step function 
can be verified by a contour integration and the theorem of residues. 

The t- integration of expression (23) can be carried out with the help of equation (24): 
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We used here the identity 


6(x) 


-i/ 

2n 


e' lx y dy 


(26) 


Substitution of equation (25) back into equation (23) gives for the retarded Green’s 
function 


GJE) = — 



AB 


(<o)(e K " /kT - 1) 


E 

co + le 


(27) 


By the same procedure the advanced Green's function is 


G„(E) 
a 2 n 



I AB ( W )(e Rw / kT - 1) 

dc o 

E 

co - ic 

R 


(28) 


These are the spectral representations of the Green's functions. 

Now consider E to be complex. Inspection of equations (27) and (28) shows that we 
can consider G„(E) and G„(E) as one function which has the properties of the retarded 
Green's function when E is in the upper half- plane and the properties of the advanced 
Green's function when E is in the lower half- plane. That is, 


G(E) 



(G r (E> 



9m E > 0 
9m E < 0 


(29) 


G(E) can be considered as a single analytic function having singularities, along the real 
axis. The retarded Green's function is analytic in the upper half-planej and the ad- 
vanced Green's function is analytic in the lower half of the energy plane! This analytic- 
ity follows from a theorem (ref. 33) which states that the complex function G(E) has 
analytic continuation in the upper (lower) half of the complex E plane if its Fourier 
transform G(t) vanishes at t < 0 (t > 0). So the functions 0(t - t') and 0 (t ' - t), which 
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act to cut off the retarded and advanced Green's function when t < t’ and t > t', re- 
spectively, are the necessary and sufficient conditions for the analytic continuations in 
the complex E plane. 

It is interesting to see the information that can be gained in crossing the cut along 
the real axis. That is, by taking the difference in the values of G(E) just above and 
below the cut, we can evaluate the spectral density, which in turn allows us to calculate 
the correlation functions. From equations (27) and (28) 


/too 

G(E + ie) - G(E - ie) I dcoI AB (a;)(e Rw/kT - 

2 71 I 

l/- co 




(30) 


Noting that the Dirac delta function can be written 


6(x) = lim 

e~0 2 tt \x + ie 



(31) 


equation (30) becomes 

G(E + ie) - G(E - ie) = iI AB ^ (eE/kT ‘ 1} (32) 

Equation (32) is an important relation. If the equations of motion of G(E) can be 
solved, we can use equation (32) to obtain the spectral density I ab (E/K), which can be 
used in equations (19) and (20) to produce 


<A(t)B(0)> =i 



G(fio> + ie) - G(hct> - ie) 
e Ko>/kT _ 1 


e RcVkT 


-io>t 

e 


du> 


<B(0)A(t)) 



G(fiqj + ie) - G(fio> - ie) e ~iwt ^ 
e ho>/kT _ 1 


(33) 


(34) 


We have, finally, in equations (33) and (34) the relations by which we can calculate 
time correlation functions from the Fourier transforms of the Green's functions. In the 
case of the simpler correlation function ( BA) , equation (34) is 
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( BA) =i 



G(Kct> +-ie) - G(Ba) - ie) dw 
e RwAT _ 1 


(35) 


The decision of how to define the Green's function for a particular physical problem will 
depend in large measure on these last three equations. One wishes the correlation 
function (BA) to be a meaningful, important property of the system under study; and 
this will dictate the choice of operators B and A, from which the Green's function 
((A(t);B(t'))> is defined. 


DISORDERED HEISENBERG FERROMAGNETS 
Equation of Motion 

Ferromagnetism is a cooperative phenomenon, and it is this fact that makes a fer- 
romagnet so interesting physically and statistically. The direction of the spin of a mag- 
netic moment is influenced by the field of other spins of the system, and the field from 
that spin in turn affects the other spins . As discussed in the section EARLIER WORK, 
the fundamental interaction responsible for ferromagnetism is the exchange interaction 
of the form J(f, g)S f • S . In this expression, S* and S are the spin operators asso- 
ciated with sites f and g, and J(f, g) is the exchange integral between these sites. For 
a ferromagnetic interaction, J(f, g) is a positive quantity. The Heisenberg model is a 
collection of magnetic moments, every pair interacting by exchange. The Hamiltonian 
for the Heisenberg model is the sum of these interactions. 

J(f,&)S f . S g -Hgju B ^S^ (36) 

2 f g f 

Since each pair of sites is counted twice in the double sum, the factor 1/2 is present. 

The second term of the Hamiltonian is the Zeeman term with the external magnetic field 
H assumed to be parallel to the z-axis, g is the Lande g-factor, and p B the Bohr 
magneton. Representing the spin operators in terms of their components 

S* - S* ± iS^ (37) 

puts the Hamiltonian into the form 
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(38) 



Z Z J(f ’ g) [l (Sf+S 8 + ^ + SfZS 3 ' H8Mfi 

f g 



s 


z 

f 


These components obey the commutation relations 


[sj,s^ = 2*8*5^ 
[ S ?,S|]= iK S f ‘6^ 


(39) 


For a system of spin- 1/2 particles where each spin is in one of two possible config- 
urations, S z = K/2 or S z = -R/2, a very suitable Green's function is one of the type 
used by Montgomery, Krugler, and Stubbs (ref. 32), ((S/(t);S“(t'))) . From this choice 

m m D 

of Green's function the correlation function (S^S^) can be calculated by equation (35). 
Since, from elementary quantum mechanics, 


S"S 


+ 


= S(S + 1)K 2 - 



- KS 


(40) 


for the spin- 1/2 case, 


S - S + 



(41) 


That is, (SjSj) is directly related to the magnetization. For spins higher than 1/2, 

there is ambiguity in the term (s z ) of equation (40) which makes ((Sj(t);S”(t'))) an 
unsuitable Green's function for spins greater than 1/2. 

The Green's function we have chosen to develop the theory of disordered ferromag- 
nets of any spin value are those used by Callen (ref. 19) in his extension of the Green's 
function theory of crystalline ferromagnets to general spin: 


Gfg(t) = 




(42) 


The t' associated with S z and S” has been set equal to zero. The superscript a is 
a parameter which will become useful only when generalizing the theory to systems of 
spin S greater than 1/2. The special case a = 0 will be the most important one when 
using the Green's function to calculate properties of the system. 


13 



To study the time development of this function in the Heisenberg model, we write its 
equation of motion 


ihA G a g (t) = h6(t)5 fg 



+ 



The averaged value of the commutator we label 0(a). 




= ©(a) 


(43) 


(44) 


Note, that for a = 0, 0 is related to the magnetization cr = (S z )/S. 


0 ( 0 ) = 



- 2R 


S 


z 

f 


Commuting (t) with the Hamiltonian allows equation (43) to be written 


(45) 


iR ~ G fg (t) " H ^B RG fg (t) + K 2 J(f,h)G f,hg (t) - RX) J(f ’ h)G h,fg (t) =K6(t)6 fg 0(a) 
^ h h 


(46) 


where 





(47) 


The appearance of higher order Green's functions Gj ^ g (t) in the equation of motion 
for G a g (t) is an expected development after the discussion of the general theory. We 
could, of course, write equations of motion for these higher Green's functions and these 
would introduce still higher order functions. Before talking about ways of decoupling 
this chain of equations we introduce the Fourier transforms of the functions in equa- 
tion (46): 
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G f, g (t) = J Gf g (E)e' i(E/K)t d ^ 




G ?,hg <‘>= /" G* hg (E)e- i(E /»‘d(|' 
J - O0 ' ^ 




( 48 ) 


5(t) = -±- 
2 




j 


Substitution of equations (48) into equation (46) allows us to express the equation of mo- 
tion in terms of Fourier components of the Green's function: 


EGfg(E) = HgM B RG^ g (E) - K ^ J(f, h)[c* hg (E) - G^ fg (E)] + - 5 - 6 fg 6(a) (49) 

h ’ 2n 

The obvious way to decouple the hierarchy of equations for the Green's functions is 
to find a way of representing higher functions in terms of lower order functions. This 
procedure presents a finite set of equations which can generally be handled. In most 
instances (refs. 16 to 21) the set of equations is cut off after the first equation by a de- 
coupling approximation that represents a second Green’s function (i.e. , one with three 
indices) in terms of the first, doubly indexed, functions. One gets the impression that 
some of the decoupling approximations which have been used were not developed a priori 
from considerations of the physics of the system but that the first criterion was to pro- 
duce an approximation that would make the mathematics tractable. Plausibility argu- 
ments can sometimes be offered as to why the approximation might be expected to be 
valid under certain conditions (of temperature, e. g. ) and the results of these calcula- 
tions do indeed compare favorably with exact results in many circumstances. But, the 
statement (ref. 18) of a decade ago that "the philosophy and justification of the decou- 
pling procedure is still far from being well understood" is only slightly less true today. 

We decouple in the manner first employed by Tyablikov (refs. 16 and 17) by making 
the approximation 


G f,hg (E) ^ < sZ > G hg (E) forf * h (50) 

which is equivalent to the statement 
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for f * h 


( 51 ) 



“ < s ?> 


aS 

(Sj(t);e g Sl 


This approximation replaces the average of a product by the product of averages. For 

aS z 

example, if f =£ g, then S? will commute with the operator e g S and 

D 




_ aS z 

\ , 

aS z 

S z (t)Sj(t),e g S' 

) = -ie(t)(s z (t) 



(52) 


In this case, equation (51) replaces the average of the product of the operators S?(t) and 
aS z - 

Sj(t),e S Sg 


by the product of the averages of these operators. This is sometimes re- 


ferred to as ignoring the fluctuations in S z , but this is not literally true. What the 
Tyablikov type of decoupling ignores is any spin correlations between sites f and h 
and between sites f and g, while uncorrelated fluctuations in S z are treated correctly. 

The effect of equation (50) is to reduce to one the types of Green's functions in the 
equation of motion 


E - Hg/igfi - K(S Z > ^J(f,h) 


Gjg(E) + K^S 2 ) 5^ J(f, h)G* g (E) = A 6 fg 0(a) (53 j 

■ 2 T[ 


Having purged the equation of motion of higher Green's functions, we can simplify equa- 
tion (53), conceptually at least, by considering it a matrix equation 


AG - 1 


(54) 


where 


r 


A fh “ 6 fh ] 


2n 


- - HgM R 
h _ 

©(a) 


(S z ) 

0(a) 


£ J(f,j) 


2tt(S z ) 

©(a) 


J(f, h) 


(55) 


and 1 is the unit matrix. 

The indices in these equations are lattice site labels; and so for a perfect crystal 
the property of translational invariance would be exploited to solve equation (53), which 
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we can do quickly. If we call the Green's function for a perfect crystal r , we would 
write its spatial Fourier transform as 




r 5f(E)=— X ™ 4 e^ f_g ^^rl(E) 
51 N Z — i k 


J(f 


,g)= iZ ei(fs)ir/(f) 


y 




5 _Z\ «i(f-g)-k 

fg 


(56) 


where k is a reciprocal lattice vector. When substituted into equation (53), these re- 
lations (56) yield 


r a (E) - 1 

k 2 tt E - E 

k 


(57) 


where 


E_ 

k 


= Hg/i B R + fi<S z )[/(0)-/(k) 


(58) 


and 


/(k)^e i(f - g) - k J(f,g) 
f-g 


(59) 


The Green's function r^.(E) has poles on the real axis at E = E__. These are the 

k k 

eigenvalues of the Hamiltonian, the energies of the elementary excitations of the system. 

From equation (32) the spectral density I^g(w) will have 5-function singularities if the 

Green's function has poles only on the real axis. In this case, from equation (19) the 

correlation function will oscillate only at the frequencies E—/R. In our case the corre- 

k 
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lation function for a = 0 is ( S S^(t)) and its frequency of oscillation is that of the corre- 

o o 

sponding spin wave. 

In a disordered ferromagnet, translational symmetry is absent and so the technique 
of expressing things in reciprocal lattice space, or k-space, is not available to us. The 
solution for the disordered case which follows uses the method of reference 32 and is 

n 

expressed in terms of T (E). 

k 


Introduction of Disorder 

Before proceeding, it would be beneficial to define more precisely what is meant by 
the term disorder in our context. We have stated that a Fourier transform into k-space 
would be disallowed for equation (53) if translational symmetry were absent. To be 
more complete, we should say that even in a symmetric lattice the Fourier transform 

n 

that yielded F __(E) in equation (57) would not be valid if the exchange interactions were 
k — — 

site dependent, that is, if J(f, g) depended on both indices and not just on r^ - r . So, 

another type of disorder which can be discussed is a randomness in the value of ex- 
change interactions. Certainly, it is difficult to conceive of a system with positional 
disorder that did not also have disorder in the J's since the exchange integral is very 
sensitive to the interatomic distance, but the converse is quite conceivable. Substitu- 
tional impurities in a perfect- crystal lattice would preserve positional symmetry while 
producing disorder in the interactions. The disorder to which we will be directing our 
attention is a randomness in the value of the exchange integrals with or without position- 
al disorder. The requirement on the position of the spin sites is that their ensemble- 
averaged positions (the ensemble average performed over systems with similar disor- 
der) possess lattice symmetry. This requirement implies that the members of the 
ensemble must be topologically equivalent. Topological equivalence means not only 
that a one-to-one correspondence of spin sites be preserved throughout the ensemble, 
but also that the number of near neighbors of each spin of each ensemble member re- 
main the same. Specifically, this study investigates the cases where the ensemble- 
averaged spin sites have simple cubic (sc), body-centered cubic (bcc), or face-centered 
cubic (fee) symmetries. 

What does it mean to take ensemble averages over systems with similar disorder ? 
A way to depict the disorder in a system graphically is shown in figure 1, where the 
number of exchange interactions per unit energy of strength J is plotted against J. 

The case of a perfectly ordered ferromagnet with nearest- neighbor interactions of 
strength J® would be represented by the delta function spike, (1/2)ZN6(J - J®), where 
Z is the number of nearest neighbors, or coordination number, and N is the number 
of spins in the crystal. The dashed curve might represent a disordered ferromagnet 
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Strength, J 


Figure 1. - Distribution of exchange interactions. 


whose exchange interactions are distributed randomly about a mean interaction, say J^. 
Those systems whose exchange distributions match this curve would be said to have 
similar disorder, and an ensemble whose members have identical curves is the kind of 
ensemble we refer to when performing our averages. We assume, also, that the disor- 
der is distributed throughout the ferromagnet and that deviations are not correlated over 
finite distances. 

It is not necessary to know the exact distribution function, rj{ J) of figure 1, to be 
able to calculate quantitatively the effects of disorder in a magnetic system . This work 
parameterizes the distribution of J-values by two quantities, a mean value of the ex- 
change integral J^(f, g) and the mean-square deviation from this mean. We shall char- 
acterize the amount of disorder in a system by a disorder parameter p which is related 
to the ratio of the mean- square deviation and the mean: 



With respect to figure 1, p might be thought of as a measure of the width-height ratio of 
the distribution function curves. Having defined disorder, the task now is the solution of 
the equation of motion for a disordered ferromagnet (eq. (53)). 

Equation (54) is the equation of motion reduced to matrix form. For the case of a 
perfect crystal this equation is 
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a °f = i 


(61) 


where the superscript on the linear operator A indicates perfect order. We define a 
matrix 


A = A - A 


(62) 


Simple manipulation of equations (54), (61), and (62) leaves 

A°G - AG = 1 (63) 

- 0 

Since r commutes with its inverse A , operating from the left with r changes equa- 
tion (63) to a Dyson equation 


G = f + f AG 


(64) 


or 


g = r + TAr + TArAr + rATArAr + 


(65) 


We wish to obtain an ensemble-averaged Green's function (Gj g (E)> , which is done by 

performing an ensemble average over systems with similar disorder. Since r^(E) is 

~ ^ k 

the Green's function for the perfect- crystal (or zero disorder) case, (T) = T; and so 
the ensemble average of equation (64) is 


( G) = f + f <AG> 


( 66 ) 


which we can put in a Dyson form 


< G> = f + f S< G) 


(67) 


where S, the "self- energy" is 


S=<^G)<G> _1 (68) 

This allows us to express ( G) in terms of f and the self- energy: 

<G> - (1 - fs) -1 f 


20 


(69) 



The quantities ( G) and 2 are translationally invariant because they are ensemble- 
averaged quantities and, as such, are diagonal in k- space. The Fourier space transfor- 
mation that effects this diagonalization yields 



cL 

Evaluation of 2__(E) will allow us the solution to the ensemble- averaged Green’s function 

a k 

<G a (E)>_ 

Upon ensemble averaging, equation (65) becomes 

(G) =r + f(A)r + r< at a) r + r < at at a> f + . . . (7i) 

where A involves the differences in the values of the exchange integrals between the 
perfect and disordered systems. If we assume the deviations of the J f g's from the 
mean to be symmetric, all ensemble-averaged quantities in equation (71) involving odd 
powers of A will vanish. 

<A> .0 

(AT Af A) = 0 
(Af Af Af AT A) = 0 


This assumption defines r as the Green's function for a perfect crystal whose exchange 
interactions are the mean interactions of the disordered ferromagnet, which we call the 
corresponding perfect crystal. 

Now, we make the approximation 

2 2° H (Af A) (73) 

A t\ __ A n n 

which in equation (71) is equivalent to replacing (A* 211 ) by (A^) . We have carried out 
the ensemble averaging of equation (73) in appendix B. In the difference matrix A, we 
designate the deviations from the mean exchange interaction as 
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( 74 ) 


j(f,g) = J(f,g)- J°(f,g) 
and the ensemble averaging involves use of the relation 


(j(f,g)j(h,i)> = j 2 ( f >g)[ 5 fh 5 gi + 5 fi 6 g h] (75) 

2 

Here j (f, g) is the mean square deviation for the exchange interaction between sites f 
and g and is proportional to the disorder parameter p (eq . (60)). Using relation (75) 
in equation (73) gives 



2r V 


+2j 2 (g, j)(r g . - r gg )> (76) 


The diagonalization of has also been performed in appendix B. In reciprocal 
lattice space the self-energy becomes 



r^,[2/ 2 (0) - 2/ 2 (k) - 2/200 + / 2 (k + k') + / 2 (k - k’)] 
k 

(77) 


As shown in appendix B, when applied to the three cubic systems with nearest- neighbor 
interactions, equation (77) simplifies to 


si(E) = 2irp 

k R0(a) 



(78) 


We can now write the Green's function for the disordered system explicitly in terms 
of the eigen energies of the corresponding perfect crystal. Using equation (78) in equa- 
tion (70) gives 
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1 


E_ 

k 


( 79 ) 


( G a (E)> _ =MW 
k 277 


E - E_ - pE_ — 

k k N E - E_, 

k' 



This is the Green's function which we will now employ to study the effects of disorder on 
a ferromagnet. 


Density of States 

Having finally solved the Green's function for a disordered ferromagnet in equa- 
tion (79), we can immediately use it to calculate the density of spin wave states g(E). 
This is a quantity of central importance in the calculations of the magnetic properties 
because it allows sums over k-states, triple integrals, to be replaced by a single inte- 
gration over the energy of these states. For example, we will later wish to calculate 
the magnetization' <j, which involves summing the quantity |exp(E—/kT) -ll”* over all 

k 's. This is generally done by an integration over the first Brillouin zone, 



where v is the volume of the unit cell. Alternatively, this sum can be performed with 
a significant saving of time and effort by 


g(E) 



dE 


(81) 



23 



The density of states here is the fraction of the total number of spin wave states per 
unit energy that have energy E . 

For a given crystal structure the energy of a spin wave in zero field depends on k 
and the temperature T. The T dependence comes in by virtue of the presence of (S 2 ) 
in the expression for the energy, and (S z ) decreases with increasing temperature. 
Dividing the spin wave energy by (S z ) produces a temperature- independent "energy" 
that depends only on k . There are computational advantages to normalizing the spin 
wave energies in this fashion since we shall be dealing with systems of various geome- 
tries and spin values. We define a dimensionless energy 


x_ 

k 



E iT 1 


h(S z ) j' 


> 


X 


E 

fi( S Z ) J' 


(82) 


For a specific symmetry, g(x) will be unique, whereas a g(E)- against- E curve would 

be valid. only for a particular S and T. J' has the same units as an exchange integral, 

_2 

(Energy) x (Angular momentum) , and is dependent on the system under consideration. 
In systems with nearest- neighbor interactions only, for example, J' has the following 
values for simple cubic, body- centered cubic, and face- centered cubic geometries. 

From appendix A, 


J sc^ 12J "^ 


J bcc = 16J 


J fcc - 16J ^ 


(83) 


where J is the nearest- neighbor exchange interaction. 

The Green's function (G^(E)_) is a natural function from which to derive g(x) since 

k 

it is also independent of temperature. To show this, we can substitute equations (45), 
(58), and (82) into equation (79). This gives for the zero-field Green's function, 
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K 


1 


x_, 

k 


( 84 ) 




<G°(x)>_ 

k 


a function which depends only on k and x for a given system and not on T or S. 

To derive g(x) formally, we go back to the Green's function of a specific disordered 
system and recast the diagonalized form of equation (54) as 


,-L W =G°(x) 


G° a (E) = ~ aB ' 

a P 2 n E R 2nJ’ R a P 

9 B aa 9R x B of 

/ nZ\ Zn 


2n ( S 2 ) 


where are the elements of the matrix 


fh 2fi 


5f h E J ( f ’ h > 


(85) 


( 86 ) 


after diagonalization. Calling the discontinuity in G at the real axis in the energy plane 
the "imaginary part" of G, that is, 

9m G(E ) = lim i[b(E + ie) - G(E - ie)l (87) 

e— 0 L J 


we note that 

- e o/3 a = f a a(J 5 (x - ^ B a J (88) 

The trace of the imaginary part of the Green's function, then, is the density of spin wave 
states 

— - Tr 9m G° (x) =- V 6 (x - ^ B U g(x) (89) 

2R N ap N \ J' “/ 

a 

Obviously, g(x) is normalized, since there are N eigenvalues B^. 
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( 90 ) 


/ BWdx ^Z/ 6 ( x ‘p B ») dx=i 

a 

When the perfect- crystal Green's function (eq. (57)) is substituted in equation (89), 
the density of states for the ordered system becomes 

g ( x ) =^I Tr 9m r°(E) - Vsfx - x A (91) 

u 2HN k V k/ 

k 

Equation (91) is evaluated for lattices of sc, bcc, and fee symmetries in appendix A. 

The evaluation involves replacing the sum over k states by a threefold integration over 
the first Brillouin zone. 

Finding the density of states for a disordered system g (x) is more complicated 

0 o p 

since ( G (E)) _ is more complicated than T (E). Using equation (79) in equation (89) 
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The self-energy term that appears in the disordered Green's function adds complexity to 
the problem in that g p (x) cannot be written as a sum of delta functions. Before summing 
the two terms inside the parentheses of equation (92), it is advisable to rewrite the sum 
in the self-energy term in a different form. Using the perfect crystal gp(x) just calcu- 
lated, 

N / x ± le - x_, 
k 

From the theory of functions of a complex variable (ref. 34), we know that 

x'g 0 (x') 
x - x r ± ie 

where the P before the integral sign indicates the Cauchy principle value. Now, equa- 
tion (92) can be written 



where 
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! ( x ’ x ic) S ^ x iT xg o< x ) 


Replacing the sum over the first Brillouin zone in equation (95) by an energy integration 
involving the density of states puts g^(x) into a form more convenient for computational 
purposes: 


pxg 0 (x) f 


x'g 0 (x') 
R(x,x')^ + I(x, x')2 


dx’ for 0 < x ^ 1 


gp(x) = { 


(97) 


J g 0 (x')6[R(x, x')] dx' for 1 < x 

where the vector sum in R(x,x') has also been replaced by an integral: 


R(x, x’) = x - x’ - px 


, p f x,, gQ (x,,) 

J X - X " 


dx' 


(98) 


Equation (97) is the density of spin wave states for a disordered ferromagnet, the 
quantity that will allow us to calculate the effect of disorder on magnetic and thermody- 
namic properties. But g (x) is important and interesting in its own right since it shows 
the manner in which the energy states are redistributed when disorder is introduced in 
the lattice. The physical insights gained from studying density- of- states curves are 
helpful in determining when disorder will and will not be important. 

Although we discuss in detail the density- of- states curves of the three cubic struc- 
tures with and without disorder in the section APPLICATION TO SPECIFIC SYSTEMS, 
there are general remarks and observations which can be made at this stage. For ex- 
ample, the shape of the density- of- states curve at low energies is important in the cal- 
culation of low-temperature properties since in Boson systems it is only these low- 
energy states that are occupied when T 0. We can get information concerning the 


28 



low-energy states in the disordered system by investigating (G®(x)) k in the limit of 
small x. Since 



the Green's function for small x can be written 


<G°(x))_ « — (100) 

k x— 0 nJ' x - x_(l - p) 
k 

The poles occur when x = (1 - p)x k . Another way of stating this is that the introduction 
of disorder has scaled the energy of the lowest eigenstates downward by the factor 1 - p. 
There are more spin waves, then, at lower energies in the disordered ferromagnet, and 
we would expect g p (x) to be accordingly larger than gg(x) for x near zero. We can be 
more quantitative about this effect without becoming too complicated. 

At low energies the dispersion relation for any crystal must be quadratic 

x_ ^ C ct k^ for small k (101) 

k Sl 


where C gt is some constant which depends on the lattice structure. In this low-energy 
domain we can calculate the density of states by integrating over a small sphere of 
radius k j : 


g 0 (x) 




1 


x— 0 4 
3 




( 102 ) 


The quadratic dispersion relation results in square- root behavior for the density of 
states at low energies. If these lowest states had their energies scaled downward by the 
factor 1 - p, 

_2 

x£ = U - p)C gt k (103) 
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what would be the effect on the density of states ? Going through the procedures of equa- 
tion (102) again yields 


g (x) « (1 - p)" 3//2 g 0 (x) 
p x-0 u 


(104) 


The effect of disorder on the low-energy density of states is to increase the density by 
(1 - p)“ 3 / 2 . This is an important result of the theory, and this factor will reappear 
when discussing low-temperature properties such as magnetization and specific heat. 
This increased number of low-energy states is a reflection of the decrease in the ex- 
change coupling between some pairs of spins that comes about with the introduction of 
disorder and the corresponding decrease in the spin wave energy. 

There will be pairs of spins in a disordered ferromagnet that are coupled more 
strongly than any in the corresponding perfect crystal. We would expect, therefore, 
that there are spin wave states in the disordered system with energies higher than any 
in the ordered system. Inspection of equation (97) does indeed show that g ^(x) can be 
finite above x = 1; whereas, gg(x) = 0 for x > 1. 

We can make additional general remarks concerning the effects of disorder on the 

XL 

density of states in regard to the moments of g (x), where the n moment of g (x) 

h' r 

is defined 



(105) 


Obviously, the zeroth moment is unaffected by p since gp(x) is always normalized by 
the way the density of states is defined. The same is true of the first moment; that is, 
the average spin wave energy does not change with disorder. We can show this with the 
help of equation (89) 


x g(x) dx = J x i- 


6(x - — B 
J’ 


Of 


dx 


2h j_ 

J’ N Z_J 


B 


a 


(106) 


a 


The average spin wave energy is seen to depend on the trace of matrix B defined in 
equation (86), a quantity depending only on the average exchange coupling which is inde- 
pendent of the disorder parameter p. 

The negative first moment, however, does depend on the amount of disorder. We 
will show later that this quantity is proportional to the reciprocal of the Curie tempera- 
ture. In the same manner as above and using equation (85), 
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— 277 Tr G°(0) 
2fi 


( 107 ) 


f I^dx = — — V J_ 
J X 2K N ZL/ B q; 

a 


The trace of the Green's function at zero energy, then, is proportional to the negative 
first moment. If we evaluate this for the ensemble- averaged Green's function (eq. (84)), 
we get 



(108) 


That is, the minus first moment varies with the amount of disorder as (1 - p)~*. 

To summarize, then, even before investigating the density of states of specific sys- 
tems and their relation to disorder, we have been able to deduce several general fea- 
tures: 

(1) We expect the low-energy portion of the g- (x)-against-x curve to get larger 

-3/2 ^ 

with increasing disorder as (1 - p) ' . 

(2) We expect g (x) to extend above the energy band of the perfect crystal. 

H 

With these changes of shape: 

(3) The area under the curve will remain unchanged. 

(4) The average energy will remain unchanged. 

(5) The negative first moment of the curve will vary with disorder as (1 - p ) ~ . 

It is interesting to see how these observations are verified later in the detailed density- 
of- state calculations. 


Correlation Functions for Spin Greater than 1/2 


Having found the Green's function for a disordered ferromagnet, the central problem 
becomes the derivation of the corresponding correlation function from which we hope to 
calculate the magnetization, or (S z ) . The correlation function that equation (34) allows 
us to calculate from the Green's function (79) is the spatial Fourier transform of 
^aS z 

{e g S g S+(t)}, which we label ^/(k,a): 
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Substitution of equation (79) into this equation would allow the straightforward calculation 
of ^/(k,a), but this correlation function is useful for the S = 1/2 case only. For ex- 
ample, for a = 0 


1^4/(k,0)e i(f - g) '^ = <S“S+(t)> 
k 


( 110 ) 


Setting t = 0 and f = g gives, for the spin- 1/2 system, the magnetization as discussed 
in equations (40) and (41): 


R < s f> = y " <s f "(o)sJ-(o)> 


(ill) 


For systems with S > 1/2, however, the magnetization cannot be calculated from 

i//(k,a) as it stands. The added complexity inherent in higher spin systems results from 

z 2 

the ambiguity in the ( (S ) ) term in equation (40), which we repeat here for convenience: 


<Sf S J> - S(S + 1)K 2 


(( s i ) 2 )-*< s f> 


( 112 ) 


This term is not a variable when S = 1/2 since it has the same value for the two pos- 
sible spin configurations. Callen (ref. 19) has developed a technique for crystalline sys- 
tems that enables one to relate (Sjs£) and 0(a) to (S^) for any spin value. We can use 
some of his results to extend the theory of disordered ferromagnetic systems to higher 
spins. The technique involves utilizing the dependence of i£/(k,a)and 0(a) on the pa- 
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rameter a and the solution of an auxiliary differential equation in a. Only the essen- 
tials of the Callen method are presented here. 

The correlation function that is of main interest is the Fourier transform of equa- 
tion (109), specifically the case where the site labels are the same and t = 0, 

/ aS ? 

/e S £ Sj 



(113) 


The dependence of the right side of equation (113) on a is contained in 0(a), which can 
be factored out to leave 



= 6(a)<f> 


(114) 


where 



(115) 


The technique is to represent equation (114) as a differential equation by writing both 
/ aS f - +\ / aS f ^ 

(e S £ S £ ) and 0(a) in terms of derivatives of fe ] with respect to a. This differen- 
tial equation is 
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2 /aS z \ fi(l + $)e aK +$ 


(1 + $)e aR - 4> 


I aS? \ 9 /aS z \ 

e M - S(S + \W (e M = 0 


(116) 


which has for its solution 


aSj\ $ 2S+l e -aRS _ (1 + $) 2S+l e aK(S*l) 
/ [$ 2S+1 - (l + $) 2S+1 ] [(1 + $)e aR - 4> 


(117) 


Finally, (S 2 ) comes from equation (117) by differentiation. 


<S Z > =-! (e aS ' 
da \ 


_ R (S - $)(1 + <i>) 2S+1 + (S + 1 + $)$ 2S+1 

n ^2S+1 -2S+1 

(1 + 4 >) - $ 


(118) 


By this equation one can solve for (S z ) self- consistently for any spin. 

The <£'s in equation (118) are the analogous quantities in a disordered system to 
Callen's perfect-crystal 4?'s, which we call We define 4> R as the sum in k-space 

of the occupation numbers |exp(E_/kT) - 1 for spin waves with reciprocal lattice 

vectors k . This can be expressed in terms of the density of states as 


. 0 1 
$ = _ 


exp U‘ 1 


/ 


g Q (x) dx 

hj’(s z )x 


(119) 


where gg(x) is the crystal density of states. Equation (115) defining $ for a disordered 
system can easily be put into this form with the help of equation (92): 


gp(x) dx 
^HJ'(S Z ) x 


= *(P, T) 


( 120 ) 
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APPLICATION TO SPECIFIC SYSTEMS 


In the preceding main section the general theory for disordered ferromagnets was 
developed. This included the derivation of a Green's function and the resulting density 
of states for disordered systems. These quantities allowed the self-consistent evalua- 
tion of the magnetization. The theory has been developed to a stage where we can now 
study in more detail the effect of disorder on specific systems. Which properties and 
which structures are affected most by disorder and in what temperature range ? Is the 
spin value important and what does the theory predict in the case of substitutional dis- 
order as in the mixed chalcogenides of europium ? These and other questions are 
treated in this section. 


Cubic Systems 

Density of states . - From equation (97) the density of spin wave states for a disor- 
dered ferromagnet can be calculated provided the density of states of the corresponding 
perfect crystal is known. We have calculated the density of states for sc, bcc, and fee 
crystals, g gc (x, 0), g bcc (x, 0)and g fcc (x, 0), respectively, and shall summarize these 
calculations here. A more detailed version of these calculations is presented in appen- 
dix A. The technique of Bowers and Rosenstock (ref. 35) was employed in these calcu- 
lations in which the density of states is written as 

g(x) =^- 6^x - xjj (121) 

k 

The dispersion relation, the energy's dependence on the wave vector k , is substituted 
in the argument of the delta function; and the sum over k- states is replaced by three in- 
tegrations in the usual way. The dispersion relations for the normalized, dimensionless 
energy defined in equation (82) are, for nearest- neighbor interactions, 
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x_ = - - - (cos ak, + cos ak + cos ak_) for sc crystals 
k 2 6 x y z 


x_ = — - — (cos - k cos - k cos - kJ for bcc crystals 

k 2 2 \ 2 x 2 y 2 




> ( 122 ) 


x_ = — - - (cos - k cos - k + cos - k cos - k„ + cos - k cos - k 
k 4 4 \ 2 X 2 y 


2 y 2 z 2 z 2 x 


for fee crystals^ 


where a is the edge length of the unit cube. 

Performing the first of the integrations implied in equation (121) results in an inte- 
grand for the second integration that can be put into a standard form of a complete ellip- 
tic integral (ref. 36) of the first kind. The final integration was performed numerically 
on an IBM 7044-7094 using Gaussian quadrature and, where functions were smoother, 
Simpson's method. 

The densities of states for perfect- crystal sc, bcc, and fee ferromagnets are shown 
in figure 2. Notice that singularities appear in all three spectra. For both the bcc and 
fee cases the density of states has logarithmic singularities, and there are discontinu- 
ities of slope in the sc and fee spectra. The slope discontinuities at x = 1/3 and 
x = 2/3 on the sc curve and at x = 3/4 on the fee curve all occur where the density of 
states has square-root behavior with slope approaching infinity. These are the familiar 
van Hove (ref. 37) singularities that appear because of the periodicity of the lattice. 

When disorder is introduced in a ferromagnet, the density of states undergoes a 
more or less continuous change of shape, so that, at small values of the disorder pa- 
rameter, gp(x) is very similar in appearance to gQ(x). Exceptions to this continuous 
change of gp(x) with p occur where there are infinities in gg(x). These infinities be- 
come finite when p becomes nonzero. The reasons for this will be discussed shortly. 
The most noticeable general effect of disorder on the density of states is to enlarge the 
number of spin wave states at lower energies, to decrease the number in the upper por- 
tion of the perfect- crystal energy band, and then to extend the maximum energy to high- 
er values. These features can be observed in the energy spectra of disordered systems, 
figure 3, where g gc (x,p), gb cc ( x >P) and gf cc ( x >P ) are plotted with respect to x for 
various values of the disorder parameter p. All the shape changes that depend on dis- 
order are caused by the self- energy term, from which all p- dependence originates. A 
closer look at this term would be useful. 

We recall that p came into the disordered Green's function (79) through a term 
which mixed the energy variable x and the eigenvalue x^, the so-called self-energy, 
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Figure 3. - Density of spin wave states for disordered ferromagnets as function of reduced energy. 
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With its linear dependence on the disorder parameter, the self-energy vanishes with 
disorder, reducing the disordered Green's function to the perfect- crystal function. Re- 
placing the sum by an integration allows us to rewrite the self- energy as 




The function A(x) is the sum over all the perfect- crystal eigenstates of the reciprocal 
of the difference in energy of the variable x and the eigenenergy x^. We would expect, 
then, A(x) to be negative for x near the low end of the energy band; to be positive for 
x near the high end of the band; and to fall off to zero like x * as x became larger 
than x = 1, the perfect- crystal energy maximum. Figure 4 shows these features of 
A(x) for the cubic lattices. The function -1 + A(x) , which contains all the x-dependence 
of the self-energy, is also plotted in the same figure. 

From these plots it can be seen that A(x) has discontinuities in slope at the same 
energy values where slope discontinuities appear in gg(x). Notice also that A(x) and 
gg(x) have infinities at the same energy value. By writing the expression for the disor- 
dered density of states (eq. (97)) in a way that mentions A(x) explicitly, it is easier to 
see the effect of the critical points of the self- energy on g n (x). 
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( 126 ) 


In the region 0 < x < 1, g p (x) appears as gg(x) modulated by an integral function 
of x. Since slope discontinuities will appear in the integrand only at energies where 
they appear in gg(x), gp(x) will exhibit discontinuities in slope at the same values of x 
but at an altered value for the state density. In this same region, the perfect- crystal 
energy band, a more dramatic change occurs where gg(x) has infinities, that is, at 
x = 1/2 in the bcc system and at x = 1 in the fee system. Since A(x) exhibits the same 
type of singularities at these points, inspection of equation (126) shows that the inte- 
grand, with the square of A(x) and g Q (x) in the denominator, becomes vanishingly 
small. This results in a depression of g^(x) near these values of x at finite values of 
p. That is, the introduction of the smallest amount of disorder causes the density of 
states near these critical points to change from logarithmically infinite behavior to a 
decay to zero. Figures 3(b) and (c) show that the energy range over which this depres- 
sion in gp(x) occurs is small for minor disorder and grows with p. 

We have observed that the maximum allowable spin wave energy becomes larger as 

disorder is introduced. The behavior of g (x) at energies above the perfect- crystal 

P _ 1 

band (i.e. , at x > 1) is related to the factor {l + p[-l + xA(x)]} , which modulates 

both the function gg and the argument of gg. There are two patterns of behavior for 
gp(x) for x > 1, one for the sc and bcc systems and the other for the fee system. In 
the former case, A(x) is finite at x = 1 and falls monotonically towards zero for in- 
creasing x such that the aforementioned modulating factor goes from a positive finite 
value to unity in the same range. The overall effect on the sc and bcc state density is 
to extend the upper portion of the spectra preserving the shape, which is of square-root 
nature. In this energy range the effect of disorder is to produce a high-energy tail, 
with the spectra showing a slight discontinuity in slope at x = 1 . 

The behavior of the state density in this higher energy region is not the same for 
the fee system. Because A(x) falls off from an infinite value at x = 1, the factor 
{l + p[- 1 + xA(x)]}~l is zero at this value of x and climbs asymptotically to unity. 

This has two effects. First, Sf cc (x, p) is depressed immediately above x = 1 because 
of this factor's appearance as a coefficient of gg. Secondly, the argument of gg in the 
right side of equation (126) is swept through the entire energy range as x goes from 1 
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to the new energy maximum. The result for the fee case is that not just the high-energy 

end of the crystal spectrum is reporduced in a modulated way but the entire crystal- 

state density reappears in the region 1 < x < x^ . although rescaled in a nonlinear 

m dx 

fashion. That is to say, in the energy region above the crystal band, gj cc ( x > P) begins 
at zero, exhibits a slope discontinuity, and finally has a logarithmic infinity at the top 
of the disordered energy band. 

An important part of the density of states is the low-energy end. Since spin waves 
act as bosons, at low temperatures only the low-energy spin wave states are occupied. 

In calculating magnetic properties at low temperatures, then, only the shape of the 
density- of- states curve at low energies is important. As has been mentioned in the 
previous chapter, disorder increases the state density by the factor (1 - p )~ 3 ^ 2 while 
preserving the square-root behavior. Since the limiting behavior of the perfect- crystal 
state densities has been evaluated in appendix A, we are able to write the low-energy 
part of the disordered density of states: 


g sc (x ’ p ) ~ ^ 

x-0 (1 . p) Z/2 2 


V* 


'A 


gbcc(x,p) x-o — 

x u (i - ppV 


g fcc (x ’ p > “ 


V* 


X ~° (1 - p) 3 / 2 7T 2 J 


(127) 


Magnetization and the Curie temperature . - The cooperative nature of ferromagnet- 
ism becomes very evident when we study the spontaneous magnetization of a ferromagnet. 
The dependence of the relative magnetization on temperature for a Heisenberg ferromag- 
net is shown in figure 5, where we have plotted ( S 2 ) /S = o for a simple cubic crystal of 
spin 1/2. The behavior near the Curie point T^ is of particular interest. Just below 
T^, the magnetization drops steeply and, in fact, the slope of the curve becomes nega- 
tively infinite at Tq. A small increase in the amount of thermal fluctuations in this 
region produces a relatively large decrease in the magnetization to zerb, an abrupt 
change from a state where the spins act in consort to the uncooperative paramagnetic 
phase. This precipitates an intriguing question. Since the spontaneous magnetization is 
very sensitive to thermal fluctuations near the Curie temperature, would the introduc- 
tion of disorder have a significant effect on the magnetic behavior in this region? The 
questions of whether the Curie temperature would be altered by disorder and whether the 
slope discontinuity in the magnetization curve at T^ remains are among those answer- 
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Temperature, T 

Figure 5. - Relative magnetization as function of temperature 
for simple cubic (sc) systems of spin 1/2. 


ed in this section, which treats the effects of disorder on magnetization over the entire 
temperature range. 

In the preceding main section the expression for the averaged z- component of spin 
was given for Heisenberg systems of any spin value. For the sake of convenience we 
shall rewrite it in the form of the relative magnetization, 


CT _ <S Z ) _ 1 (S - $)(1 + $) 2S+1 + (S + 1 + $)$ 2S+1 
S ~S (1 + ^S+l _ *2S + 1 

where $ in turn is a function of T, S z , and p, 



(128) 


(129) 


so that equation (128) is a self-consistent expression for the magnetization. 

Upon expansion, equation (128) has the form of the quotient of two power series, 


a =■ 


1 + a^$ + + • • • + a 2S-l^^~* 

1 + b +-b 2 $ 2 + . . . + b 2g <f> 2S 


(130) 
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where the constants a p (i = 1,2, . . . , 2S - 1), and bj, (j = 1,2, . . . , 2S), depend on 
the value of S. Before presenting the results of equation (130) for systems of various 
spin, structure, and disorder, we are interested in the behavior of a for two limiting 
cases: T approaching zero and T approaching T^, from below (i. e. , T-T^-). For 
low temperatures it will be important to know the coefficients of the lowest power of 4>, 
aj. and bp and for T near the coefficients of the higher powers of $ will be 
needed. These are 




a l =' 


2S + S - 1 


bj = 2S + 1 


J 


(131) 


and 


a 2g _i =- (2S 2 + 3S + 1) 
3 




b 2g = 2S + 1 


a 2S _ 2 =- (4S 3 + 4S 2 - S - i) b 2g _ 1 = S(2S + 1) 
6 


a 2S-3 = ~ (4s4 ' 5s2 + ^ b 2S-2 S ( 2S + 1)( 2S - !) 


(132) 


Low temperatures: When the temperature approaches absolute zero we expect all 
ferromagnets, no matter how weak or strong the exchange interactions, to have all their 
spins alined parallel. In the absence of thermal fluctuations there will be no spin re- 
versals, and the ground state of a disordered ferromagnet will exhibit the identical 
magnetization as the corresponding perfect crystal. If, in the manner of Dyson (refs. 11 
and 12) we were to expand a in a power series in T, the leading term would be the 
same for both ordered and disordered ferromagnets, namely unity. The question which 
is of interest is where in this series does the effect of disorder first appear and how 
does the disorder parameter enter. It will be shown that all subsequent terms are af- 
fected by disorder. 

At low temperatures, <i> is a small quantity because of the rapid exponential decay 
of the integrand as the energy increases. The factor jexp (E_/kT) - lj” 1 in the inte- 
grand of $ has an infinite, integrable peak at the zero of energy which, at low temper- 
atures, serves to weight the density- of- states function heavily at the low-energy end. 
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Because the low- temperature magnetization is insensitive to the shape of the high-energy 
portion of g(x), we are justified in using the low-energy form of g ^(x) in the evaluation 
4> in this region. Thus, for H = 0, 



(133) 


where C le is a constant whose value for various structures can be found from the low- 
energy forms of the state density in equations (127). The right side of equation (133) can 
be evaluated by recognizing that it can be put into an integral form of the Riemann zeta- 
function £(s) by a change of variable: 


S(s) =• 


’ <s) i * y -i 


dy 


(134) 


This gives for the low-temperature form of 4> 


$ C=r 

T— 0 


7T 1/2 C le ?(3/2) / \3/% 


kT 


2(1 - p) 3 ^ 2 \hJ'(S z )y 


That is, for T near zero, $ goes as T 3//2 and (1 - p) 3//2 . 
Equation (130) in the limit of low T becomes 


(135) 


a ^ 1 + (a, - b, )<f> = 1 - — <f> 
T— 0 1 1 S 


(136) 


Substitution of equation (135) allows a to be written to a first approximation as 


t » 1/2 c le ;(3/2) / kT ^ /2 

T-0 2S 5 / 2 (1 - p) 3 / 2 \h 2 JV 


(137) 
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which is Bloch's well-known T*^ 2 law revised to include the effects of disorder. High- 
er terms will also have dependence on p, but near T = 0 the effect of disorder is to 
increase the T 2 ^ 2 term by (1 - p)“ 2//2 , which means that the magnetization of a disor- 
dered ferromagnet will decrease more quickly than the corresponding perfect crystal in 
this temperature range. 

Curie temperature: The energy of a spin wave is not uniquely determined by its re- 
ciprocal lattice vector but depends also on the magnetization. As the temperature ap- 
proaches the Curie point from below and the magnetization drops to zero, there is a 
corresponding decrease in the spin wave energy. To evaluate $ in this limit, we shall 
expand the exponential term of the integrand in a power series. If we set the constant 
quantities of the exponent equal to W, 


w - ft 2 s J ’ 

k 


(138) 


then 




(139) 


Because the quantity Wox/T is small near the Curie temperature, the term in the 
brackets in equation (139) can be rewritten as a power series that gives 

♦ - JL f 

T-T c - Wo lx 2 


f / 


xgp(x) dx + 


(140) 
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We shall write an expression for 1 /ct, the inverse of equation (130); and since $ 
is large near the Curie temperature, only its higher powers will be kept: 


1 ~ 
a T-T c - 


, *2S , .2S-1 L . .2S-2 

b 2S^ + b 2S-l $ b 2S-2^ 

0 .2S-1 0 i2S-2 i 2S-3 

a 2S-l $ + a 2S-2^ + a 2S-3^ 


(141) 


The quotient on the right side of equation (141) is a power series in ct/T containing only 
odd integer powers: 


I « 

a T-T c - 




(142) 


More than the leading terms were kept in both the numerator and the denominator of 
equation (141) in order to verify the absence of even powers of ct/T and in order to be 
able to evaluate the constant in equation (142), which is 


B _ (2S - 1)(2S + 3)W 
1 20(S + 1) 



(2S 2 - 2S - 1)W 
4(2S + 1)(S + 1) 



xg p (x) dx 


(143) 


The Curie temperature can be evaluated from equation (142). Rearranging terms 
gives 


1 


1 - 


CT 


3T 

(S + 1)W 



dx 




(144) 


Note that this equation has no solution for arbitrarily large T. For finite cr, the right 
side of the equation remains finite as T — 00 , while the left side does not. It is an equa- 
tion valid only to the temperature where ct vanishes. As ct, and therefore the right 
side of the equation, approaches zero' the left side will become infinite unless the ex- 
pression in the brackets becomes zero. This condition determines T^: 
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dx 


(145) 


To put kT^ in terms of the exchange interaction energy, 


kT C _ S(S + 1) 
h 2 j’ 3 



gp(x) 

X 


1-1 


dx 


(146) 


As mentioned in the previous main section the Curie temperature is inversely pro- 
portional to the negative first moment of the density of states. And it was shown in 
equations (107) and (108) that 



(147) 


This implies that a disordered ferromagnet has a Curie temperature that is smaller 
than the T^ for the corresponding perfect crystal by the factor 1 - p. 


( T c)p ■ (I - p, ( T c) p= „ 


(148) 


The Curie temperature, then, is linear with the mean-square deviation of the exchange 
integrals . 

The calculation of T^, is not complete until the integral of equation (27) is evalu- 
ated. Machine calculations of this integral have been performed, using the three state 
densities derived earlier. Our results agree with the work of Watson (ref. 38), who 
calculated similar integrals analytically. Thus, we have 
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/ 

/ 


g sc (x, 0 ) 


dx = 3.033 


8bcc< x -°> 


dx = 2. 786 


dx = 1.793 


( 149 ) 


Finally, substituting for the structure-dependent factor J’ by means of equation (83), 
and using equations (147) and (149), the Curie temperatures for disordered sc, bcc, and 
fee structures of any spin value can be written in units of the nearest- neighbor exchange 
energy: 


4S(S + 1)(1 - _ p) _ j 3 jg g(g + i)(i _ p) f or sc structures 
3.033 



< 


16S(S + 1)(1 - p) 
3(2.786) 


= 1.914 S(S + 1)(1 - p) 


for bcc structures (150) 


16S(S + 1)(1 ~ P ) = 2 . 975 S(S + 1)(1 - p) for fee structures 
L 3(1.793) 


Although disorder has a significant effect on the Curie temperature, the shapes of 
the magnetization- against- temperature curves will be similar near T^,. From equa- 
tion (144) it can be seen that the magnetization of a disordered ferromagnet will have the 
same square-root behavior as the perfect crystal with a negatively infinite slope at T^. 
To first order, equation (144) can be written 


a ^ 
T-T, 



(151) 


Intermediate temperatures: Having investigated the spontaneous magnetization of 
disordered ferromagnets at low temperatures and near the Curie point, we shall present 
in this section the results of calculations of <j over the entire temperature range of the 
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ferromagnetic phase. These calculations have been carried out for various values of 
spin and disorder and for the three cubic structures. 

The method used to generate the various a- against- T curves is as follows. The $ 
corresponding to some ct/T ratio was evaluated by a machine integration of equa- 
tion (129). Equation (130) was then used to solve for a for any desired spin system. 

To find the temperature corresponding to the value of ct, we divide by the original ct/T 
ratio. To show the effect of spin in the theory, these calculations were carried out for 
spins of 1/2, 1, and 7/2. The S = 7/2 case was included for two reasons. First, it 
provides an example of behavior for a system of spin significantly higher than the com- 
monly studied S-values. Second, the europium chalcogenides, probably the best exam- 
ples of Heisenberg ferromagnets, have spin values of 7/2 associated with the europium 
ions; and it is desirable to have calculations that are applicable to real systems. The 
appropriate equations for these values of S are 



Figure 6 shows the magnetization curves for three perfect- crystal ferromagnets of 



0 1 2 3 4 5 6 

Temperature, T, f/j/k 


Figure 6. - Average z-component of spin (s z ) as function 
of temperature for simple cubic (sc), body-centered 
cubic (bcc), and face-centered cubic (fed systems for 
spins of 1/2 and 1. 
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sc, bcc, and fee symmetries. Two values of spin, S = 1/2 and S = 1, are presented 
for each structure. In this figure the average z- component of spin is the ordinate rather 
than the relative magnetization a. For the same value of spin the Curie temperature is 
seen to be highest for the fee crystal, intermediate for the bcc, and lowest for the sc. 
This might be expected, since a spin in a fee lattice interacts more strongly with the 
rest of the system, being involved with 12 nearest neighbors. The magnetic moments 
associated with ions in bcc and sc lattices interact with eight and six neighbors, respec- 
tively, and have correspondingly lower Curie temperatures. 

The effect of the value of S can also be seen in figure 6. For a given structure the 
Curie temperature increases with S as S(S + 1). Since the exchange interaction be- 
tween two spins (J^Sj • Sj) increases as S , this increase in T c is understandable. 

The effect of disorder on the spontaneous magnetization is shown in figure 7 for a 
sc system of S = 1/2, a bcc system of S = 1, and a fee system of S = 7/2. The be- 




Fitjure 7. - Relative magnetization as function of temperature for disordered cubic systems of various spins. 
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havior in each case is qualitatively the same. At any finite temperature, the magneti- 
zation of a disordered ferromagnet is smaller than for the less disordered system; and 
the Curie point in each case decreases linearly with 1 - p. 

The a- against- T curves for all the systems considered, independent of the structure 
or of the value of S or p, all start from a value of unity for a at T = 0 and decrease 
monotonically to zero where there is a slope discontinuity at T = T c - The behavior of 
these curves near T = 0 and T = T c is similar (although the displacement of T c due 

to disorder is significant), having T 3//2 and T 1//2 dependence, respectively, with the 
coefficients of these terms being functions of p. Any changes in shape of these curves 
due to disorder, then, are subtle changes not easily discernible from the plots as they 
have been presented. Changes in shape due to p would be more easily studied by 
plotting a against T for a system with various amounts of disorder in such a way that 
the curves intersect at T = 0 and T = T c - That is, the temperature scale would be 
altered for each value of p to allow the T c 's to be superimposed. This procedure has 
been followed in figure 8(a), where bcc systems with S = 1/2 and with various amounts 
of disorder are presented. Displayed in this fashion, it can be seen that the effect of 
disorder is to flatten the curve shape. The change in behavior of the magnetization 
curves brought on by disorder are not just changes in scaling. This flattening effect of 
disorder is less strong at higher spin values, as can be seen from figure 8(b), where the 
same type of plot is done for a bcc system of S = 7/2. Recently, Sharon and Tsuei 
(ref. 39) have confirmed this flattening phenomenon in measurements on amorphous 
ferromagnetic Fe- Pd- P alloys. 

Specific heat . - To calculate the magnetic or spin contribution to the specific heat 



(a) Spin 1/2. (b) Spin 7/2. 

Figure 8. - Relative magnetization as function of rescaled temperature for disordered body-centered cubic systems 
of spin 1 12 and 7/2. 
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CL, one needs an expression for the thermal average of the Hamiltonian. Mills (ref. 40) 
has derived such an expression exactly from spectral functions of the type discussed in 
connection with the general theory. Cooke (ref. 41), by an alternate method, has also 
derived a thermal average (3C) which although not exact for S > 1/2 is simpler in form 
than the Mills result. The higher order correlation functions that Cooke ignored, how- 
ever, do not contribute to the leading terms of (3C) at low temperatures. Thus, his re- 
sult is a good approximation for any spin in this temperature region. Cooke's expres- 
sion is 




E°(k) + - 2(k,T) 
2 



(153) 


where Eg is some temperature-independent constant that does not affect the specific 
heat, E°(k ) is the spin wave energy of a magnon at T = 0, 

E°(k) = R 2 s[j(0) - J(k)] (154) 

and Z(k , T) contains the temperature dependence 

S(k,T) = fi^S 2 ) - fis][f(0) -/(k)] (155) 

Changing the summation in equation (153) to an integration gives 
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At low temperatures, the integral in this expression is easily computed by a change 
of variable and use of equation (134) 


N 


E 0 + - 
u 2 


<S Z > + ^(S 2 ) 


C le /3 tt 1//2 \ ? /S\/_J^_\ 5/2 

(1 - p) 3//2 \ 4 / 


(157) 


For T near zero, equation (157) indicates that the leading temperature term for (3€) is 
of order T^ 2 and has dependence on disorder as (1 - p)" 3 ^ 2 . Since the specific heat 
is given by 


C 


s 


d(3Q 

dT 


(158) 


it will behave as (1 - p)'3/2r ^3/2 a ^. i ow temperatures. Disorder, then, will have im- 
portant effects on low-temperature thermodynamic properties of ferromagnets, effects 
that should be easily detected experimentally. For example, a value of 0. 1 for p pro- 
duces a 17 percent increase in the low-temperature specific heat over the value of the 
corresponding perfect crystal. 

Paramagnetic phase . - The cooperative effects of ferromagnetism disappear above 
the Curie temperature. In this region the system of permanent magnetic moments goes 
into the paramagnetic phase, where the ordering tendency of the exchange interaction is 
subdued by the competing influence of thermal fluctuations and the net magnetic moment 
is zero in the absence of external fields. To investigate how disorder affects the para- 
magnetic phase of a Heisenberg system, an expression for 1/a such as equation (141), 
is needed. For the sake of computational simplicity, the S = 1/2 case is treated and 
the general results are shown to be valid for any S. For S = 1/2, equation (141) be- 
comes 


= 1 + 2 $ 


(159) 


The right side of equation (159) can be put in the form of a hyperbolic function: 


1 + 2 
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If we write the part of spin wave energy with external field dependence as h and the 
other energy term as y(x), the right side of this equation can, by the use of a mathe- 
matical identity, be rewritten so that 



where 


h = HgMgli and 


y(x) = (K 2 /2 )ctJ'x. 


Setting 


t n = tanh — — 
u 2kT 

► 


t, = tanh 
1 2k 


and bearing in mind that only t^ has dependence on x, 
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By use of the series expansion for the hyperbolic tangent, t^ can be written as 
series in reciprocal temperature: 
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where C n is the n n moment of the state density 
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n 



dx 


Equation (164) can be used to solve for a iteratively to any order in 
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(162) 


(163) 


a power 


(164) 


(165) 
For ex- 
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ample, where a = tg is used as a zeroth approximation, the second approximation gives 



In this equation, it can be seen that the effect of disorder does not appear until the third 

-2 

term, where Cg appears in the coefficient of T . We have shown previously that Cj, 
the first moment of the density of states, is related to the average energy of spin waves 
and is unaffected by disorder. The second moment C 9 and higher moments of g (x) 
will have p dependence, but they do not appear in the leading terms of the series. 

The paramagnetic Curie temperature 9 is a parameter of experimental interest 
and it is defined by the Curie- Weiss law, 


X = 


c 

T - 6 


(167) 


where C is the Curie constant and x is the magnetic susceptibility. The paramagnetic 
Curie temperature is found empirically by extrapolating the linear high- temperature re- 
gion of a (x)~ ^-against-T curve to the (x)" 1 = 0 axis. Generally, 9 is appreciably 
greater than the actual ferromagnetic transition temperature T^,. Since the suscepti- 
bility is proportional to cr/H, we can find 9 by multiplying equation (164) by H, keep- 
ing only lower order terms in T"*: 




- (gMgR)” 1 — 


1-^1 


*0 

\4kT/ 


For high temperatures, tQ h/2kT and so 


[X ]" 1 « (g MgR )" 1 



when [x] 


= 0 , 


then 


C 1 fi 2 J’ 

T =— 

4k 


(168) 


(169) 


(170) 
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The paramagnetic Curie temperature for a system with disorder is identical to the 
perfect- crystal result. 

The conclusion to be drawn from this section is that the paramagnetic phase of a 
Heisenberg system is less sensitive to disorder than the ferromagnetic phase. We 
might have intuitively expected this, since in the paramagnetic region the thermal fluc- 
tuations have become sufficiently large to destroy the cooperative effects. With this 
much "thermal disorder" already present in the system, we might have suspected that 
the additional disorder of randomness in exchange interactions would be less significant 
than in the lower temperature phase. 


Europium Chalcogenide Mixtures 

Divalent europium can be made to react with the chalcogenide series, O, S, Se, and 
Te to form crystals of the cubic NaCl structure. In 1961 (ref. 42), EuO was discovered 
to be ferromagnetic; and soon afterwards, ferromagnetism was found in EuS and EuSe, 
while EuTe was found to be antiferromagnetic (refs. 43 and 44). These compounds are 
insulators and are nearly ideal Heisenberg ferromagnets . Only the Eu ++ ion is mag- 
netically active in these compounds, and it has a spin of 7/2. As the size of the non- 
magnetic ion and the lattice parameter increase, the ferromagnetic coupling decreases 
and so the Curie temperature decreases from 69 K for EuO to 16. 5 K for EuS to 7 K for 
EuSe. The antiferromagnetic EuTe has a Neel temperature of 7. 8 K. 

Mixtures of these chalcogenides of europium will be topologically equivalent, but 

the exchange integral between two Eu ++ ions will vary according to which chalcogenide 

ions are nearby. A system disordered in this fashion should be able to be treated by 

the theory developed in earlier sections. We shall consider the EuX^Z^ ^ system 

where X and Z represent any of the chalcogens - O, S, Se, or Te. Our convention 

will be to let tj be a measure of the concentration of the "weaker" specie, that is, 

(T^J > (T r ) . The assumption will be made that the X and Z ions are ran- 

V VeuZ V VEuX 

domly distributed in the lattice. In the NaCl structure the Eu ions form a fee sublattice 
and the nonmagnetic ions form another fee sublattice. Each Eu ion has six nonmagnetic 
nearest neighbors and, in its own sub lattice, 12 Eu "nearest neighbors. " We assume 
that exchange interactions exist only among nearest- neighbor Eu ions, and we pose the 
following problem. What is the Curie temperature for such a system, and how does 
Tq depend on the relative concentration rj ? 

To answer this question we propose a simple model. There are two nonmagnetic 
ions that lie closest to a line joining two nearest- neighbor Eu ions, and the size of these 
two ions will be most important in determining the strength of the exchange integral. 

We assume, then, that the exchange integrals will have one of three values, depending 
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on which types of ions occupy the two sites closest to the Eu-Eu bond: 


J = J : 
J = J 2 

j =!(j + j ) 

2 1 z 


when two X ions occupy these sites 
when two Z ions occupy these sites 

> 

when an X ion and a Z ion occupy these sites 


(171) 


In order to calculate the disorder parameter p, which will be a function of 77, we 
will need to know Jg, the average exchange interaction or the exchange integral of the 
corresponding perfect crystal. If N is the number of Eu ions, the total number of ex- 
change interactions is 6N and 


Number of interactions of strength Jj = 6N77 

o 

Number of interactions of strength Jg = 6 N (1 - 77) 

Number of interactions of strength - (J! + J 0 ) = 6N2t)(1 - tj) 

2 1 1 


> 


(172) 


The average exchange integral is, therefore, 

J 0 = tj 2 Jj + (1 - T7) 2 J 2 + 277(1 - 77) 1 (Jj + J 2 ) 

2 

= 77 Jj + (1 - 77 ) J 2 


( 173 ) 


The disorder parameter for this system is 


P(i 7) = 


! «J - J 0 ) 2 ) 


6 


J 


2 

0 


- 1 d 2 77(1 - 77) 
12 (dT? - l ) 2 


(174) 


where d is the ratio 
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(175) 


d = 



1 


From equation (150) the Curie temperature is proportional to the product of 1 - p and 
J Q . Therefore, the Curie temperature for EuX Z, that is T r (r»), in units of the 

T c for pure EuZ, is V 


T cfr) 

T c (°) 


= [1 - P(r?a - 

Jo 


1 - d T ?(! - v) 

12 (d r] - l) 2 


(1 - d 77) 


(176) 


This relation has been used to calculate the Curie temperatures of the three sys- 
tems - EuS tj O^_ tj j, EuSe^j ^ and EuSe^O^ ^. The results are shown in figure 9. 




Figure 9. - Curie point as (unction of relative 
concentration (or europium chalcogenide 
mixtures. 
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When the value of d is relatively small, that is, when the Curie temperatures for pure 
EuX and pure EuZ are relatively close, the change in Tq(t]) with 77 is almost linear. 
This is so because p(tj) never gets very large when d is small, in which case the 
change in T ^( 77 ) is caused almost entirely by the change in Jq(t]), which is linear. If d 
becomes sufficiently large, the [l - p(rj)] term assumes more importance, producing a 
departure from linearity. In fact, for d> 0.923 a minimum occurs in the T c (t})- 
against - 77 curve. In this case, replacing the X ions of EuX with Z ions would initially 
cause a reduction in even though EuZ is a stronger ferromagnet with a higher T c 

than EuX. For d > 0. 98 the ferromagnetic phase is destroyed completely by small con- 
centrations of Z (i.e., 7? <: 1). For the EuSe^O^ ^ system, d = 0.90, a value suffi- 
ciently large to warrant an experimental investigation of the region of low oxygen con- 
centration to see if any of these surprising predictions associated with high d- values are 
observable in this system. 


DISCUSSION AND CONCLUSIONS 

Several general statements can be made now, in answer to the central question posed 
in the INTRODUCTION. There it was asked whether disorder would have a significant 
effect on a cooperative phenomenon like ferromagnetism and if so, how? The Green's 
function theory developed herein to study distributed disorder in a Heisenberg ferromag- 
net gives us an affirmative answer, and the disorder has been shown to manifest itself 
in several ways. 

Disorder, introduced by allowing a randomness in the strength of the exchange 
couplings, has a marked effect on the density of spin wave states. Modest values of the 
disorder parameter p can produce relatively large changes in the state density g^(x) 
in the form of enhancement of the low-energy densities and extension of the energy band 
to higher values. Physically, this can be explained as an effect of the addition of great- 
er and lesser exchange interactions than are found in the corresponding ordered system . 
Some spins in the disordered system will find themselves more strongly coupled to their 
neighbors, while other spins will be less strongly coupled than in the perfect crystal. 
Since the energy of a spin wave is proportional to this exchange coupling, there will be 
more very low-energy magnons and also higher energy magnons than the case where all 
spins are equally coupled. This is analogous to the situation in lattice dynamics 
(ref. 45) where the substitution of light impurities into a crystal produces localized 
modes in the phonon density of states outside the perfect- lattice continuum, while heavy 
defects produce an enhancement in the lower part of the spectrum. Both of these effects 
are present in the spin wave state density, where allowing a randomness in the values of 
exchange interactions corresponds to adding heavy and light impurities to a phonon sys- 
tem. 
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The increase in the density of low-energy states is of the order (1 - p)~^^ , and 
this manifests itself in the low-temperature properties such as the specific heat. When 
T is near zero, all but the lowest spin wave states are unoccupied. As the disordered 
ferromagnet is heated slightly, the spin waves are allowed to occupy the states to a 
higher energy according to the Bose distribution, but there are more of these states 
than in the ordered system. Thus, more energy must be put into the system to produce 
the equilibrium distribution associated with a given temperature, and the specific heat is 
correspondingly larger. 

The spontaneous magnetization of a disordered ferromagnet decreases with temper- 
ature more quickly than for a crystal and the Curie point was shown to diminish linearly 
with disorder as 1 - p. The transition from the ferromagnetic to the paramagnetic 
phase is no less abrupt for a disordered system, however, in that the spontaneous mag- 
netization exhibits the same negatively infinite slope at T^,. 

In the preceding main section, the effects of disorder on systems of various struc- 
ture and spin were presented. Without repeating specific results, there are general 
conclusions to be drawn from this work. Probably, the most important statement to be 
made is that order and disorder are of most significance in cooperative phenomena and 
of less consequence where cooperative effects are absent. A corollary to this statement 
might be that the effects of disorder are more noticeable at low temperatures than at 
high. The disorder parameter is involved in the first term of the low- temperature 
specific- heat expansion, but it is not present in the high- temperature susceptibility until 
the third term of the series. 

At the conclusion of what has been an initial attempt to include the effects of dis- 
tributed disorder in the theory of magnetism, there are suggestions that can be made 
for future work in this area. Certainly, the decoupling approximation used in the sec- 
tion Equations of Motion limits the accuracy of our calculations. Improvement in the 
low-temperature expansions (ref. 19), for example, are sure to come about by the use 
of more sophisticated decoupling schemes that might include correlation effects. An- 
other improvement to the present work is likely to be produced by a more complete de- 
scription of the disorder. The one parameter description used herein, treats all dis- 
tributions of exchange interactions alike, provided the relative mean-square deviations 
are equal. However, the shape of such distributions are undoubtedly important for some 
properties. And, of course, to produce results that agree more closely with experi- 
mental measurements, inclusion of more than just nearest- neighbor interactions might 
well be required. 
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In conclusion, we have generalized the theory of ferromagnetism to include the ef- 
fects associated with the disappearance of translational symmetry. It is our hope that, 
in rectifying some of the incompleteness in our knowledge of magnetic systems, we have 
provided new insights in our understanding of all cooperative phenomena. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 28, 1972, 

502-01. 
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APPENDIX A 


DENSITY OF SPIN WAVE STATES OF PERFECTLY 
ORDERED CUBIC FER ROM AG NETS 

In calculating the density of spin wave states of a disordered ferromagnet, use is 
made of the density of states of the corresponding perfectly ordered ferromagnet. The 
details of the derivation of the density of states for crystals of simple cubic (sc), body- 
centered cubic (bcc), and face centered cubic (fee) symmetry are presented here. Sim- 
ilar calculations have been carried out by Jelitto (ref. 46) and others (ref. 47) whose 
results are in agreement with those presented here. 


Simple Cubic 

According to Green function calculations, a spin wave of wave vector k has an 
energy 


E(k) = Hgp.^ + fi<S z )[/(0) - /(k)] 
When there is no external field, this dispersion relation becomes 


E(k ) = R(S Z ) 


I 

f-g 


fgL 


1 - e 


-i(f -g )■ k 


(Al) 


(A2) 


where J ^ is the exchange interaction between the spins at sites f and g. The vectors 
f and g are the position vectors for these sites . Exchange interactions are relatively 
short ranged and so we make the assumption 



if f - g = Nearest-neighbor vector 
otherwise 


For a simple cubic lattice, then 


(A3) 
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(A4) 


E(k) = ^(S 2 ) J 



-iA- k 
e 


where A represents the six possible nearest- neighbor position vectors ±ax, ±ay, ±az, 
where a is the length of the cube edge and x, y, and z are unit vectors in the x, y, 
and z directions. Carrying out the summation in equation (A4) gives 

E(k) = 2h(S z ) J(3 - cos ak - cos ak - cos ak ) (A5) 

x y z 

It is convenient to deal with a "reduced energy, " which is the energy normalized to 
range between zero and unity. We define, then, 



- - - (cos ak + cos ak + cos ak ) 
2 6 x y z 


(A6) 


We wish to derive a spin wave density- of- states function g fx) which represents 

Ot 

the number of spin wave states of energy x per unit energy, divided by the total num- 
ber of states. The fraction of states with energy between Xj and Xg, then, would be 

X2 g sc (x)dx (A7) 

1 



From the way we described g (x), it is obvious that g (x) is a normalized function 

SC sc 



(A8) 


The reciprocal lattice vectors k which are associated with the spin waves of energy 

x k are uniformly distributed throughout k- space. So the derivation of g sc (x) reduces 

itself to the calculation of the volume bounded by the two surfaces of constant energy: 

x(k , k , k) = x and x(k , k , k ) = x + dx. The fraction of the total volume bounded 
X V z x y z 

by these two surfaces is g (x) dx. 

Since any spin wave can be represented by one wave vector in the first Brillouin 
zone (i. e. , any spin wave with a wave vector outside the first Brillouin zone is indis- 
tinguishable from the spin wave with the corresponding wave vector in the first zone), 
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we need only concern ourselves with the first Brillouin zone when wishing to sample all 

—► O 

k -space. The volume of the first Brillouin zone is (n/a.)' . We can write g c (x), then 

SC 

as 



(A9) 


where 6 is the Dirac delta function. 

Before calculating g (x) from equation (A9), it might be beneficial to say a few 
words concerning the appearance of the delta function. If one were interested in count- 
ing the number of spin wave states in the energy range between two energies x^ and 
Xg, one way would be to use the delta function in the following manner: 



6(x - x k ) dx = Number of states in energy range 


Xj < X < Xg 


(A10) 


As the sum is made over all possible k- vectors, the delta function acts as a counter, 

registering unity every time k is such that x 1 < x_ < x 9 and registering zero when 

1 k z 

the wave vector is such that its spin wave energy falls outside this range. Comparing 
equation (A10) with expression (A7) suggests that 


Ssc^Z 6 ^-^ (All) 

k 

The number of spin wave states is of the order of magnitude of the number of atoms in 
the crystal. In macroscopic samples, this number is sufficiently large to warrant re- 
placement of the sum over k- states with an integration over the Brillouin zone. Thus, 


N 



(A12) 


With this as the rationale for equation (A9), it remains only to carry out the three 
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integrations. After substituting equation (A6) into equation (A9) and making use of the 
identity 


6[f(x)] = 


df(x) 


dx 


6(x - x Q ) 


(A13) 


where f(xg) = 0, and making the change of variables 


cos (ak x ) = u 


cos (aky) = v 


cos (ak ) = w 


(A 14) 


we have 


g 


sc 




5[w - (3 - 6x - u - v)] 
■y/(l - u 2 )(l - v 2 )(l - w 2 ) 


(A15) 


We shall see that after performing the w- integration, the limits of integration in the 
u-v plane must be redefined. That is, since 


5(w - 3 + 6x -i- u + v) 

y (i - u 2 ki - v 2 )(i - w 2 ) 


dw =-Z 


y(l- u 2 )(l - v 2 ) 1 - (3 - 6x - u - v) 2 


for - 1 < (3 - 6x - u - v) < 1 


0 otherwise 


(A16) 


the original limits of the u and v integrations must now be more strictly defined so 
that not only -1 < u < 1 and -1 < v < 1 but also 

u + v > 2 - 6x and u + v < 4 - 6x (A17) 

After the first integration, then 


g 


sc 




^(1 - u 2 )(l - V 2 ) 1 - (3 - 6x - u - v) 2 


(A18) 
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E-6936 



where for 0 < x ^ 1/3 


for 1/3 < x < 2/3 


and 


and for 2/3 < x < 1 


H 

ii 

2 

V U = 1 

u^ = 1 - 6x 

V L = 2 “ 

Uy = 3 - 6x 

V U =1 

U L = 1 

V L “ 2 “ 

U U ~ 1 

V U = 4 " 

u^ = 3 - 6x 

V L “ ” 1 

u^ = 5 - 6x 

V U = 4 ' 

U L “ -1 

1 

li 

> 


6x - u 


6x - u 


6x - u 


6x - u 


The second integration, involving the inverse square root of a fourth-degree poly- 
nomial in v, can be reduced to an elliptic integral (ref. 36) of the first kind. After 
factoring the integrand of equation (A18), we have 


r u 

I dv 

^(v - l)(v + l)[v - (4 - 6x - u)][v - (2 - 6x - u)] 


(A19) 


The factors in the integrand of expression (A19) can be arranged in such a way to look 
like 


I 

[(v - a)(v - b)(v - c)(v - d)]' 1//2 


(A20) 


where 
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a > b > c > d 


The roots a, b, c, and d are, of course, 1, -1, 4 - 6x - u, and 2 - 6x - u arranged 
to satisfy the inequality (A20). When so arranged expression (A19) takes the form 



dv 

y (v - a)(v - b)(v - c)(v - d) 


(A21) 


Such an integration, between the third and second largest roots of a fourth-degree poly- 
nomial, can be reduced to a complete elliptic integral of the first kind 


K(m) = 



V 


dt 

(1 - t 2 )(l - mt 2 ) 


where the argument m is 


_ (a - d)(b - c) 
(a - c)(b - d) 


and the coefficient of K(m) is 


(A22) 


(A23) 


A — (A24) 

y<a - c)(b - d) 

Carrying out the second integration in this manner leaves 



K(m) 

|/ (1 - u 2 ) 


du 


(A25) 


where 




The integrations involved in equation (A25) were carried out on an IBM 7094-7044 
and the results are shown in figure 2(a). As could have been predicted by the symmetry 
of earlier equations, the density of states is symmetric about x = 1/2. There are two 
discontinuities in the slope of the curve, occurring at x = 1/3 and x = 2/3. The be- 
havior immediately below x = 1/3 and above x = 2/3 is of square-root nature, which 
means that the slopes have infinities of inverse square-root nature. By well-known 
arguments the low-energy density of states must be proportional to the square root of 
the energy and the proportionality constant can be calculated by evaluating equation (A25) 
in the limit of x approaching zero: 



It is convenient to have an analytic expression of the density of states in the low-energy 
limit because, in the calculation of low- temperature thermodynamic quantities, the low- 
energy portion of the density-of- states curve is the only part that has importance. 
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Body-Centered Cubic 


The density of spin wave states for a body- centered cubic ferromagnet g^ cc (x) is 
calculated in this section in a manner analogous to the procedure for the simple cubic 
case. The zero-field dispersion relation is written by carrying out the sum indicated 
in equation (A2), again considering only nearest- neighbor interactions. After defining 
a reduced energy, the delta- function formalism introduced for the simple cubic case is 
used to find the volume in k- space between surfaces of constant energy, a process in- 
volving a threefold integration over k- space. A plot of the results of these integrations 
is presented in the form of g^, cc (x) against x, and some comments are made on the 
shape of this curve. 

In the body- centered cubic system, each spin has eight nearest neighbors. The dis- 


persion relation resulting from summing terms like 
bors is 


J, e 

fg 


i(f -g)-k 


over nearest neigh- 


E (k ) = 8K< S 2 ) J [1 - cos — k cos - k cos - k \ (A27) 

\ 2 x 2 y 2 z / 

where, as before, J is the exchange integral between nearest neighbors and a is the 
length of the edge of a body- centered cube. The maximum energy that a spin wave in a 
bcc system can possess is, then, 16h( S z ) J. Dividing equation (A27) by [E(k)]„ 0 „ 
yields a reduced energy that is convenient in that it is a dimensionless quantity which 
varies between zero and unity and has the added advantage of removing the temperature 
dependence (( S 2 ) is temperature dependent) from the spin wave energy. For the bcc 
system, then, the reduced energy is 

x_ = - - - cos - k cos - k cos - k (A28) 

k 2 2 2 x 2 y 2 z 


Setting 


cos - k. = u 


cos - k = v 
2 y 


cos - k = w 
2 z 


(A29) 


we have 
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(A30) 


SbccW'-T 



dw ■ 


5 w - 


1 - 2x 


uv 


j uv 


|/(1 - u 2 )(l - v 2 )(l - w 2 ) 


The w- integration will produce nonzero results only when u, v, and x are such that 


-1 < \ . Z 2 - g < l 
uv 


(A31) 


Inequality (A31) together with our knowledge that u and v can only have values between 
+ 1 and -1 fixes the limits of the u and v integrations after the w- integration is com- 
pleted. The first integration leaves 


Sbcc w - 4 



dv — ■ ±. - ■ 

|/(1 - u 2 )(l - v 2 ) uV - (1 - 2x) 2 ] 


(A32) 


Since the integrand is an even function of u and of v, it is sufficient to integrate in one 
quadrant of the u-v plane only because of the symmetry of the integration limits. 

Following the procedure of the simple cubic case, the v- integration can be put into 
the form 


Sbcc (x > = 



dv 


^(v - a)(v - b)(v - c)(v - d) 


(A3 3) 


where for 0 < x < - 
2 


o1 , 1 - 2x „ l-2x H n 

a = 1 b = c d = - 1 


u 


u 


and for - < x ^ 1 
2 
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= 1 b = -J— lx c = j_- 2x d= .j 


u 


u 


which gives 


g bcc^ ~ ~ 



K 


du 


' u - 1 1 - 2x j 
lu + 1 -2x 


(u + 1 1 - 2x ) 




(A3 4) 


The density of states generated by equation (A34) is shown in figure 2(b). The in- 
teresting portions of the curve, near x = 0 and x = 1/2, have been calculated analyti- 
cally; and the other points have been evaluated by numerical machine integration. As in 
the simple cubic case, the curve is symmetric about x = 1/2. However, whereas, 
g sc (x) was everywhere finite; g bcc (x) has an infinity at x = 1/2 of squared logarithmic 
nature. We can show this by evaluating expression (A34) in the limit of x — 1/2. We 
define a parameter f which is two times the displacement of x from the singularity: 

f = 1 - 2x (A3 5) 

As x approaches 1/2 from below, f approaches zero from above. For x < 1/2, then 


8bcc (x > =1 ! 



du 


K 

If" - fl 



L\u + f/ J 



(u + f)' 




(A36) 


When f is small, the argument of the complete elliptic integral K is very close to 1 
in the range of integration, except at the very low end where u is near f. Since K(m) 
can be expanded in a power series in m^ = 1 - m, 


K(m) = ag + a^m^ + agm ^ + . . . + ^>q + b 


jin^ + i>2 m l + 


. .) lnJ- 

/ m. 


when m is near unity, the approximation 
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(A3 8) 


K(m) a; a Q + bg In — 
m l 

becomes a good one. The approximation used in calculating the limiting form of g bcc (x) 
near the singularity is 




(A3 9) 


where 


a Q - 1.38629 



Carrying out these integrations leaves 


g bcc (x) f “ 0 c o + c l ln ’ f+c 2< ln f)2 


(A40) 


where 



(a 



0 - b Q In 4) In 2 = 0.247926 


(a Q - b Q In 2) = 0.536521 



= 0. 129006 
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The relative error introduced by the approximations of equation (A39) can be shown 
to decrease as f approaches zero because the absolute error remains bounded as x 
approaches the singularity although g^cc^O does increase. 

Finally, the low-energy form of g^^x) is evaluated by noticing that the argument 
of the elliptic integral in equation (A34) is close to zero in the range of integration when 
x is small. We write, then, 


g bcc (x) “ — 
DCC x— 0 ^3 



(A41) 


Face- Centered Cubic 

The dispersion relation for spin waves in a fee ferromagnet with nearest- neighbor 
interactions is 


E(k ) = 4R( S 2 ) J(3 - cos - k cos - k - cos - k cos - k - cos - k cos - k ^ 

\ 2 x 2 y 2 y 2 z 2 z 2 x / 

(A42) 

where the symbols are defined earlier in this appendix. Dividing equation (A42) by the 
maximum spin wave energy gives for the reduced energy 

x_ - — - — 1 1 cos - k cos - k + cos - k cos - k + cos — k cos - k ] (A43) 

k 4 4 \ 2 x 2 y 2 y 2 z 2 Z 2 x ) 


Following the procedures outlined earlier, 


%cc (x) = 




dk 5 
z 





dw ■ 


5 w - 


3 - 4x - uv 


u + v 


A 

uv 


|/(1 - U 


2 )(1 - v 2 )(l - w 2 ) 


(A44) 
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After carrying out the integration over w, the v- integration is accomplished after 
putting the resulting integrand into a standard form for the complete elliptic integral of 
the first kind, leaving 


*W X) 


%cc< x > 


where 


8fcc <x) =-T I *■ 


K(m,) 


L-2x 


^(1 - u 2 )(3 - 4x + u 2 ) 


for 0 < x < 1/2 


n- (1 - x) 



2x-l 


du 


K(m 2 ) 


-2x 


V 


1 2 
1 - u 



ft I K(m 1 ) 

+ — I du 1 for 1/2 < x < 3/4 

■ ^(1 - u 2 ' /0 -- 2 ' 


y (A45) 


)(3 - 4x + u ) 


*\4x-3 

16 f K ( m 3 } 


»2x-l 


V4x^ 


(1 - u 2 )[(4x - 2) 2 - 4u 2 ] tt 3 (1 - x) 


du 


K(m 9 ) 






1 - u" 


du 


K(m,) 


V (1 - u 2 )(3 - 4x + u 2 ) 


for 3/4 s x < 1 


_ u 2 - (2x - l) 2 
m., = - — 

1 9 

+ 3 - 4x 
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m -u 2 + (2x - l ) 2 
4(1 - x) 2 


m 3 = 


4(1 - x/ 


-u 2 + (2x - l) 2 


Figure 2(c) shows the results of machine calculations of equation (A45). Unlike the 
sc or bcc cases, the density of spin wave states for the fee system shows no symmetry 
about x = 1/2. The most striking difference in the fee curve is the logrithmic singular- 
ity in g(x) at the energy maximum . At x = 3/4 there is a discontinuity in the slope 
reminiscent of the behavior of g„_(x) at x = 1/3 and x = 2/3. That is, the slope of 
gfcc(x) becomes infinite as x approaches 3/4 from below, and the behavior of 
is of square- root nature in this region. The limiting form of gf cc (x) as x approaches 
zero is easy to calculate after noticing that m ^ is small throughout the range of inte- 
gration when x is small. Thus, 



g fpc (x) « — Jx 

lcc x-0 n 2 V 


(A46) 
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APPENDIX B 


DERIVATION AND APPLICATION OF SELF-ENERGY TERM 
Derivation and Diagonalization of the Self-Energy 

The self-energy matrix is defined by equation (73), 


S a = ( AT A) (Bl) 

where A is the difference between the A matrices for the disordered ferromagnet and 
the corresponding perfect crystal, 


A 



(B2) 


and the angular brackets indicate an ensemble average over systems with similar dis- 
order . 

The components of S a in lattice space can be written 


r 


where, from equation (55), 



s !r(£2>gh r M i i i ) 
h i 




(B3) 


(B4) 


Here j(g, f) is the deviation of the exchange coupling between sites g and f from the 
mean value: 


j(g,f) = J(g,f) - J°(g,f) 


(B5) 


The superscript zero indicates that the term is associated with the corresponding per- 
fect crystal. Substitution of equation (B4) into equation (B3) gives 



2 

2?7< S Z ) -KEE j(g> h)j(j,i)(r g . 

0(a) /hi 


hi 




(B6) 
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We assume deviations of the exchange interactions to be symmetric about the mean 
value. The ensemble averaging, then, is carried out by use of the relation 

(j(g,h)j(i,j)) = j 2 (g,h)(6 g .6 h . + 6 gj 6 M ) (B7) 

2 

where j (g, h) is the mean- square deviation of the exchange integral between spins at 
sites g and h. Since j (g, h) is an ensemble-averaged quantity, it is translationally 
invariant and depends only on the ensemble- averaged |r - | . Equation (B7) assumes 

that there are no three-or-more-site correlations in the deviations. For example, the 
ensemble average of the product j (i, j)j(i,h) vanishes. Thus, equation (B6) becomes 



Next, we write equation (B8) in a more symmetric form and employ the Fourier trans- 
form of r gh , 


r 




r e Hg-h)-k 
k 


(B9) 


where k is a reciprocal lattice vector. 




j 2 (g,h)- 26 gj ^] j 2 (g,h)e i(g - h) ‘ k 
h 




(BIO) 


It will be helpful in deriving the diagonalized form of 2 to study the product of 2 
and r . 


(fg >fg=E r fg S Bi 

g 


(Bll) 
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Since 



k 


k 


(B12) 


Equation (Bll) becomes 


iff| r( M )yyy r -' r - ei 

lg \NQ(a) / 4 / j/ j k k 


(f-g)-k 


g k' k 


X ] 26 gj X j 2(g > h) _ 25 gj 2 j 2 (g> h )e l(g " h) 'k + j 2 (g,j) 


,i(g-j)-k 


+ e i (3-g)' k 


- 2 j 2 (g, h)[ 


(B13) 


An example of how this expression will be treated can be given by looking, say, at the 
second last term. Forgetting the coefficient for the moment, we have 

yyy r_,rj 2 (g,i)e i <^«>'V (f -s^' = yy r _ ir _ e i(tH)f , y J 2 (g;j)e i(i-g)(k + k-) 
< << k k I k k V 


g k’ k 


k ' k 


(B14) 


2 

Since j (g, j) is translationally invariant and depends only on (g - j), its Fourier trans- 
form is 


/ 2 (k) = ^ j Z fe>i) e ' 

g 


, 2 <<r no-i(g-j)-k 


(B15) 


which allows equation (B14) to be written 


ZZ r r r p* (f+fV<H) ' r 


(B16) 


k' k 


79 



Similar manipulation of the remaining terms simplifies equation (B13) to 


(fs) fg = 


' 2r(S z ) \ 

V N0(a)/ 

k ' 



T 

r^,r^e i(f 'i )-k [2/ 2 (0)- 2/ 2 (k) 


2/ 2 (k’) + + k ') + / 2 < k - k ' ')] ( B1? ) 


Since the matrix (f S) is diagonalized by the Fourier transform 


(ff) fi = -V(f 


NZ— / k 

-*• 1 

k 


(B18) 


comparison with expression (B17) indicates that 


/ z \ 2 

(f £) k ' = (l^) N Z r iT ' r k [ 2/2<0) ' V2<f ,) ' + /2,ir+ ?> + /i** - *'>] 

k (B19) 

The f and 2 are both diagonal in reciprocal lattice space, so that 

(fsz, = r^, (B20) 

k k k 

The self- energy can now be written in diagonal form, 


/ z \ 2 

NZ r k [ 2/2<0) ' 2/2<f '* " 2/2<ir> + /2< f + * ’> + / 2 f - * 1 ’>] 


(B21) 


Application to Cubic Systems with Nearest -Neighbor Interactions 


If we rewrite the sum in equation (B21) as 
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(B22) 


i S r iT { 2 l/2 (0) ' A* ’>] + 2 l/2 (0) ' >] 

k 

- [/2< 0 > - / 2 <* + *'>] - U<°> - - * ’>]} 

we recognize terms similar in form to the zero- field spin wave energies. For example, 
from equation (58), when H = 0 

E^=K<S z )[ / /(0)- < /(k)] (B23) 

From the definitions of ^OO and </&), equations (B15) and (59), respectively, we see 
that when only nearest- neighbor interactions are considered 

.2 j 2 E- 

S 2 (0) -> 2 (k) = j- [/(0) -y(k)] = (B24) 

J Jl^S 2 ) 

2 

Here j is the mean- square deviation for nearest- neighbor exchange interactions, and 

J is the nearest- neighbor exchange integral. We can use the results of appendix A, 

where E_ is evaluated for the cubic systems, to evaluate equation (B24). 
k 

For the face- centered cubic system, for example, 
y 2 (0) - ,/ 2 (k) = 4j 2 sin 2 ^ (k x + k y ) + sin 2 £ (k x - k y ) + sin 2 ^ (k y + k z ) 

+ sin 2 ^ (k y - k z ) + sin 2 ^ (k z + 1^) sin 2 ^ (k z - k x ) (B25) 

where an alternative expression to equation (A 42) has been employed to represent E_ 

k 

in terms of the squares of sine functions. This permits equation (B21) to be written 
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where 


'N 


k l (k x ±k y ) 


k 9 = - ( k „ ± k J y 


4 ' y z 


k 3 = \ < k z ± k x^ 


The sum in equation (B26) can be rewritten 



Because of the equivalence of the lattice sites in a fee system 


(B27) 


(B28) 



allowing expression (B28) to be simplified to 



(B29) 


(B30) 
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But 



(B31) 


again because of symmetry considerations. This further simplifies expression (B28) to 



(B32) 


which, when substituted back into equation (B26), gives 


2 * = 
k 



x_, 

k 



(B33) 


Since the disorder parameter is 


P 



(B34) 


where Z is the number of nearest neighbors 


s a 2 np E 

k h©(a) k N 



k 


(B35) 
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The other two cubic systems, sc and bcc, can be handled in the same manner since 
the symmetry arguments are the same . Doing this verifies that equation (B35) is a gen- 
eral result for the three cubic systems. 
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